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13.  ABSTRACT 


The  objective  of  this  research  program  is  to  develop  computer  programs, 
using  the  finite  element  method,  to  predict  stresses  and  deformations  in  the  vicin¬ 
ity  of  underground  excavations.  The  computer  programs  should  allow  for  arbitrar 
initial  stresses  in  rock,  arbitrary  shape  and  size  of  the  opening,  any  given  se¬ 
quence  of  construction,  non-homogeneous  material  properties,  interaction  of  rock 
with  supporting  structures,  progressive  damage  and  time  dependent  deformation 
and  load  development  on  supporting  structure.  Limited  experimental  work  to  veri 
fy  key  points  in  the  theory  is  planned. 

Research  during  the  first  year  is  directed  towards  literature  survey,  seleC' 
tion  of  mathematical  models  for  behavior  of  rock  and  development  of  computer  pro 
grams  for  elastic-plastic,  elastic-brittle  rock  and  for  progressive  failure  of  rock 
around  underground  openings. 

At  this  reporting,  selection  of  mathematical  models  has  been  completed. 
The  report  presents  finite  element  computer  programs  for  i)  Plane  Strain  Analysis 
of  Elastic- Plastic  Mohr-Coulomb  Materials  anci  iii  Two-Dimensional  Analysis  of 
a  No  Tension  System.  Plans  for  further  requirements  of  these  programs  in  line 
with  research  objectives  are  discussed. 
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Technical  Itoport  Sum mary 


I’lu*  oMivtlw  of  thl >  research  program  Is  to  develop  computer  programs, 
u  si  mu  do  If  nlti  element  method,  to  predict  stresses  and  deformations  in  the  vicin¬ 
ity  of  niuleru round  excavations,  TIu  computer  programs  will  h.ivc  the  capability 
to  ill  m  for  arbitrary  initial  stresses  in  rock,  arbitrary  shape  :md  size  of  the 
opening,  im  ulien  sequence  of  construction,  nonhomogcncous  mnterinl  properties, 
inti  ruction  of  rock  with  support  in  u  structures,  progressive  damage,  :ind  time 
dependent  deform  ill  on  :md  load  development  on  sup|>ortlng  structure.  Limited 
experimental  work  lo  verify  key  (mints  in  the  theory  is  pkinned. 

Ilesenreli  during  the  first  ye:tr  is  directed  towards  survey  of  literature 
on  the  subject,  selection  of  mathematical  models  for  mechanical  Imhavior  of  rock, 
and  development  of  computer  programs  for  clastic-plastic  Mohr-Coulomb  materials 
for  brittle  rock  following  Griffith's  theory,  ami  for  progressive  deformation  and 
fracture  of  rock  around  underground  openings  under  stress  changes  associated  with 
excavation. 

At  this  reporting,  selection  of  mathematical  models  for  elastic-plastic  Mohr 
Coulomb  materials  and  for  elastic-brittle  materials  failing  according  to  Griffith 
theory  has  been  completed.  Chapter  1  of  the  report  describes  the  theoretical  con¬ 
sider  Ulons  leading  (o  the  model  selected.  The  stress-strain  relations  for  incre¬ 
mental  or  rate  type  theory  of  plasticity  arc  generally  based  on  the  normality  rule 
and  convexity  and  regularity  of  the  yield  surface  in  a  'stress-space'.  Using  these 


concepts,  various  Investigators  have  proposed  conflicting  constitutive  equations. 

In  tills  report  tin*  clastic-plastic  behavior  of  materials  has  been  re-e.\;imlned  ns 
,i  m  itlienintlenl  generalization  of  observations  on  :i  onc-diniensionnl  test.  The 
role  of  kinematic  constraints  upon  yield  conditions  has  been  studied  and  adequacy 
of  certain  postulates  examined.  Current  theories  of  elastic- plastic  behavior  are 
found  to  be  inadequate  as  it  is,  in  general,  not  possible  to  satisfy  the  'normality' 
rule  as  well  as  ccrtinuity  of  stress  path  under  plane  strain  conditions.  Further 
research  into  this  aspect  of  material  behavior  is  needed  to  clear  the  air.  Experi¬ 
mental  phase  of  the  research  program  is  being  planned  with  this  requirement  in 
view.  In  the  mathematical  model  selected  as  the  basis  for  development  of  com¬ 
puter  programs,  a  modification  of  the  yield  surface  is  introduced  to  eliminate 
discontinuity  in  stress  paths.  For  behavior  of  elastic-brittle  rock,  the  mode! 
selected  assumes  elements  of  rock  to  be  incapable  of  supposing  tensile  and 
shearing  forces  across  a  crack.  A  review  of  literature  showed  errors  in  simil;.  r 
formulations  by  other  investigators.  These  have  been  corrected  in  the  present 
development.  The  mathematical  models  of  elastic-plastic  and  elastic-brittle  rock 
have  been  incorporated  into  finite  element  computer  programs  for  analyses  of 
stresses  and  deformations  of  plane  strain  systems.  Chapters  III  and  IV  of  the 
report  present  two  computer  codes  along  with  relevant  description,  instructions 
for  usage  and  illustrative  examples  for: 


i.  Plane  Strain  Analysis  of  Elastic- Plastic  Mohr-Coulomb  Materials 
ii.  Two-Dimensional  Analysis  of  a  Non-Tension  System 

Further  work  on  these  computer  programs  is  continuing.  However,  even  in  their 
present  form,  program  capabilities  include  consideration  of  arbitrary  initial 
stresses,  arbitrary  shape  of  openings  with  or  without  linings,  and  considerable 
variation  in  material  properties.  These  computer  programs  should  be  of  immediate 
application  to  a  variety  of  problems. 

Adequate  mathematical  models  of  rock  behavior  have  been  chosen.  The 
finite  element  method  has  been  used  successfully  to  develop  computer  codes  for 
analysis  of  complex  problems  of  stresses,  deformations  and  fracture  in  rock. 

The  method  appears  to  be  suitable  for  further  development  to  realize  the  objectives 
of  the  current  research  program. 

Experimental  work  so  far  has  been  directed  towards  development  of  suit¬ 
able  laboratory  material  (exhibiting  elastic-plastic  behavior).  No  equipment  has 
so  far  been  purchased  under  the  contract.  However,  procurement  of  a  plane  strain 
testing  machine  has  been  initiated.  It  is  expected  to  be  received  in  September  1971. 


PREFACE 


The  terrestrial  crust  is  in  a  complex  state  of  stress.  Underground 
excavations  in  this  stressed  medium  profoundly  influence  the  distribution  of 
stress  which  in  turn  determines  the  stability  of  the  opening  and  of  the  rock  in 
the  vicinity.  Traditional  methods  based  upon  the  theory  of  linear  elastic  solids 
are  inadequate.  It  is  necessary  that  the  sequence  of  construction  and  realistic 
material  properties  be  taken  into  account  in  calculation  of  stresses  and  defor¬ 
mations  in  rock. 

The  objective  of  the  present  research  program  is  development  of  finite 
element  techniques  to  predict  stresses  and  deformations  in  the  vicinity  of 
underground  excavations  allowing  for  arbitrary  initial  stresses  in  the  rock, 
arbitrary  shape  and  size  of  the  opening  .  any  given  sequence  of  construction. 
The  procedures  will  allow  for  nonhomogeneous  material  properties,  interaction 
of  rock  with  supporting  structures,  progressive  damage ,  time  dependent  defor¬ 
mation  and  load  development  on  supporting  structures.  Limited  experimental 
work  to  verify  key  points  in  the  theory  is  also  planned.  The  entire  program  is 
expected  to  extend  over  three  years. 

This  is  the  first  semi-annual  progress  report  covering  the  period  2/1/71 
to  7/31/71.  The  main  activity  in  this  period  has  been  a  review  of  the  work  done 
by  other  investigators  in  mathematical  simulation  of  mechanical  behavior  of 
elastic-plastic  solids  and  of  jointed  rock.  Models  of  stress- strain  behavior 
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have  been  proposed  and  work  in  development  of  relevant  computer  programs 
started.  A  sequential  approach  has  been  followed  whereby  a  basic  program  is 
coded  and  then  modified  to  include  all  the  ramifications  of  material  behavior  and 
actual  loading  sequences.  Two  computer  programs,  viz, 

i.  Plane  Strain  Analysis  of  Elastic- Plastic  Mohr-Coulomb  Materials 

ii.  Two-Dimensional  Analysis  of  a  No-Tension  System 

arc  included.  The  present  capabilities  of  each  program  are  indicated  in  the 
program  descriptions.  Further  development  on  all  these  is  continuing  and  will 
be  included  in  future  reports. 

The  work  is  supported  by  the  U.S.  Government  through  the  Advanced 
Research  Projects  Agency,  ARPA,  and  its  agent  the  Bureau  of  Mines,  Department 
of  the  Interior.  At  the  Ohio  State  University  the  work  is  under  direct  supervision 
of  Professors  T.H.  Wu,  R.S.  Sandhu,  and  J.R.  Hooper.  Messrs.  S.W.  Huang, 

R.D.  Singh,  C.W.  Chang  and  T.  Chang,  graduate  students  in  the  Department  of 
Civil  Engineering,  have  contributed  to  the  research  reported.  Dr.  William  Karwoski 
of  the  Spokane  Mining  Research  Center,  Spokane,  Washington  is  the  Project  Officer 
designated  by  the  sponsor.  In  early  stages  D.\  Syd  Peng  of  Twin  Cities  Mining 
Research  Jenter,  Twin  Cities,  Minnesota  acted  as  the  Project  Officer. 

The  opinions,  findings  and  conclusions  expressed  in  the  report  are  those  of 

the  authors  and  not  necessarily  those  of  the  U.S.  Bureau  of  Mines,  Department  of 

the  Interior  or  the  Advanced  Research  Project  Agency. 

R.  S.  Sandhu 
Project  Supervisor 
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CHAPTER  I 


THEORETICAL  CONSIDERATIONS 


Chapter  I.  Theoretical  Considerations 


1.1.  Mechanical  Behavior  of  Rock 

Figs.  1-1  and  1-2  show,  respectively,  typical  stress-strain  plots  for  a 
granite  and  a  marble  (Swanson  1970).  Upon  loading  the  stress-strain  curve  is 
almost  linear  and  reversible  over  a  short  portion.  Unloading  from  higher  loads 
does  not  coincide  with  initial  loading.  This  characteristic  along  with  rate  inde¬ 
pendence  distinguishes  elastic-plastic  behavior.  Reloading  closely  follows  un¬ 
loading  until  the  previous  maximum  is  reached;  whereupon  the  original  curve 
is  followed.  This  leads  to  some  simplifying  assumptions. 

i.  A  yield  point  exists  below  which  the  material  is 
linear  elastic. 

11.  The  yield  point  corresponds  to  the  maximum 
stress  level  previously  attained. 

iii.  Unloading  and  reloading  paths  are  linear,  coin¬ 
cident  and  parallel  to  the  initial  elastic  loading 
curve. 

Fig.  1-3  shows  this  simplification.  Clearly  the  yield  point  can  be  described  by 
the  permanent  or  irrecoverable  strain  or  the  area  bounded  by  the  loading  curve, 
the  unloading  curve  and  the  horizontal  axis.  Whereas  in  generalization  to  the 
three-dimensional  case,  the  stresses  and  strains  become  second  rank  tensors 
and  are,  therefore,  unordered  ,  the  area  is  still  a  scalar  product  and  retains 
Ha  ordering  characteristics.  To  this  extent,  it  is  often  preferred  as  a  measure 
of  the  elastic  limit. 
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Mechanical  behavior  of  rock  under  polyaxial  state  of  stress  has  been 
examined  In  the  light  of  brittle  failure  theories  (Brace,  1984|  Bieniawski,  1907, 
1909;  Brady,  1969,  1970).  Four  regions  of  behavior  are  identified  in  Figure  1-4. 
The  first  region  corresponds  to  closure  of  pre-existing  open  cracks  and  is  pecul¬ 
iar  to  compressional  loading.  In  region  n  material  behavior  is  linear  elastic. 
Fracture  initiation  occurs  near  the  end  of  this  region  in  accordance  with  Griffith 
or  modified  Griffith  Theory.  This  stage  also  corresponds  to  onset  of  nonlinearity 
in  the  stress  to  volumetric  strain  curve  (Brace,  1966).  Stable  fracture  propaga¬ 
tion  characterizes  region  III.  In  region  IV,  unstable  fracture  propagation  results 
in  strength  failure  and  rupture.  Differences  in  loading  and  unloading  behavior 
:ire  observed  (Walsh,  1965). 

We  have,  thus,  two  general  approaches  to  the  characterization  of  stress- 
strain  behavior  of  rock.  One  follows  the  theory  of  elastic-plastic  solids  without 
consideration  of  micro-mechanics  of  the  system.  The  other  uses  Griffith  theory 
or  modifed  Griffith  theory  to  relate  deformation  and  failure  to  initiation  and  pro¬ 
pagation  of  fracture.  It  has  been  observed  (Swanson  1970)  that  Mohr-Coulomb 
failure  lav/  applies  for  moderate  values  of  confining  pressure  and  that  at  low 
confining  pressures,  failure  is  by  rupture.  Contrary  to  plastic  behavior,  strength 
of  material  drops  to  almost  zero  in  the  direction  normal  to  the  crack  if  rupture 
theory  is  followed.  Figs.  1-6  and  1-6  depict  typical  relationships  of  failure  strength 
and  post  failure  behavior  in  relation  to  confining  pressures.  It  is  reasonable  to 
assume  that  the  material  is  linear  elastic  upto  yield  or  rupture  ,  as  the  case  may 
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Figure  1-3.  Idealization  of  Elastic- Plastic  Stress-Strain  Behavior  for  Rocks 


Figure  1-4.  Typical  Axial  and  Lateral  Stress-Strain  Behavior  of  Brittle  Rock 


Figure  1-5.  Proportional  Limit  in  Shear  for  Westerly  Gramme 


Figure  1-6.  Effect  of  Confining  Pressure  on  Stress-Strain  Behavior  of  Limestone 


be,  and  that  the  post-failure  behavior  only  is  governed  by  the  theory  used  to 
define  failure.  Both  the  elastic-plastic  Mohr-Coulomb  failure  theory  and  the 
Griffith  theory  have  been  used  in  the  course  of  present  research  to  develop 
computer  programs  for  analysis  of  stress  and  deformation  in  rock. 

1.2  Stress-Strain  Relations  for  Elastic- Plastic  Solids 

Several  approaches  have  been  used  for  formulation  of  elastic-plastic 
behavior.  Excellent  presentations  of  the  theory  are  available  in  literature 
(Drucker,  1951;  Naghdi,  1960;  Green  and  Naghdi,  1965;  Koiter,  1953;  Hill, 

1950) .  Specializations  to  Mohr-Coulomb  solids  under  plane  strain  (Drucker 
and  Prager,  1952;  Drucker,  Gibson  and  Henkel,  1955;  Reyes,  1965;  Reyes  &  Deere, 
1966)  have  been  presented.  Not  intending  to  survey  the  entire  range  of  dif¬ 
ferent  philosophies,  we  present  briefly  a  discussion  of  elastic-plastic  behavior 
under  kinematic  constraints  leading  to  the  plane  strain  formulation  for  Mohr- 
Coulomb  materials  used  in  the  computer  program  in  Chapter  HI.  We  present 
a  generalization  of  observations  on  the  c  nventional  uniaxial  test  to  the  case 
of  three-dimensions  before  discussing  the  role  of  constraints. 

An  isothermal  system  undergoing  infinitesimal  deformations  is  of 
Interest  to  the  present  report.  Extension  to  more  comprehensive  situations 
is  direct. 

1.2.1.  A  Generalization  of  One- Dimensional  Test  Results  to  Three-Dimensional 

Theory 

Analogous  to  the  case  of  one-dimensional  test,  we  assume  the  existence 
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of  a  set  D  consisting  of  all  admissible  states  of  stress.  D  is  convex,* **  has  boun¬ 
dary  B  and  includes  the  original  or  reference  state.  Clearly  D  is  six-dimensional 
(assuming  symmetric  stress  tensors,  i.e.  absence  of  body  couples)  and  is  con¬ 
tained  in  the  six-dimensional  linear  vector  space  V  spanned  by  the  six  components 
o-^j  of  the  stress  tensor.  The  stress  in  uniaxial  test  is  given  by  a  real  number  and 
in  this  case  D  is  ordered  and  convex.  To  Introduce  an  ordering  in  the  slx-dimon- 
sional  stress  space  so  that  'increase'  and  'decrease'  of  stress  are  meaningful, 
a  mapping  g  is  defined 

g:  V  — R  (1-1) 

with  the  following  properties 

i.  The  image  of  D  under  g  is  a  positive  interval 
I  C  R. 

ii.  Image  of  complement  of  D  in  V  is  the  complement 

of  I  in  R  and  g  maps  the  interior  D{  of  D  to  the  (1-2) 

interior  of  I. 

ili.  g  maps  the  'original'  or  'reference'  state  to  zero 

€  I. 

Then  D  is  ordered  by  its  image  i.e.  for  if  g(i^)  =  gj,  g(i>2)  -  82*  v2  18 

greater  than,  equivalent*«to,  or  less  than  v\  depending  upon  g£  being  greater  than. 


*Convexity  of  D  is  the  property  that  +  (1  -a)  [0,ll. 

In  literature,  there  are  frequent  references  to  a  convex  yield  surface.  This  is  in¬ 
accurate.  It  is  easy  to  see  that  convexity  of  D  does  not  imply  convexity  of  its  boun¬ 
dary  B.  Indeed,  B  is  in  general  not  convex. 

**The  term  equivalent  is  used  because  g,  in  general,  is  not  one  to  one. 

Thus,  Vj,  \>2  may  be  distinct  while  their  images  coincide. 
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equal  to  or  less  than  gj.  The  mapping  preserves  convexity  of  D.  The  interval 
I  =  [o,f]  where  f  is  the  image  of  boundary  B  of  D. 

In  one- dimensional  tests,  the  set  of  admissible  stress  states  may  include 
negative  points.  The  limiting  stress  states  in  tension  and  compression  provide 
a  positive  supremum  and  a  negative  inf imum  to  this  set.  The  boundary  B  of  this 
set  is  clearly  discontinuous.  In  a  multi-dimensional  stress  space,  B  may  be 
continuous.  To  ensure  correspondence  between  the  one- dimensional  and  multi¬ 
dimensional  cases,  g  is  a  two  to  one  mapping  in  the  case  of  one-dimensional 
loading.  Thus  the  image  of  B  is  in  all  cases  the  supremum  of  non-negative 
interval  I. 

The  boundary  B  of  D  and  hence  its  image  f  under  g  are  defined  by  prior 
deformation  and  load  history.  Considering  components  of  the  strain  tensor, 
on  the  analogy  of  the  results  of  one-dimensional  test,  a  plastic  strain  tensor 
with  components  c  is  defined  such  that 

i.  For  a  given  e"ki,  there  Is  one  to  one  correspondence 
between  elements  of  D  and  a  set  of  points  in  the  six¬ 
dimensional  space  spanned  by  e'kl  =  ekl  -  e"kl  •  e'kl 
are  identified  as  components  of  the  elastic  strain  tensor , 

d-3) 

ii.  If  a  generalization  of  Prandtl's  simplifying  assumption 
is  admitted,  the  one  to  one  correspondence  between 
<rki  E  D  and  e  is  independent  of  prior  deformation 
history,  and  for  stress  states  defined  by  interior  Dj  of 
D,  *"kl  vanish. 

A  positive  measure  of  history  of  deformation  can  be  defined  in  various 
ways.  If  plastic  strain  is  used  as  representative  of  deformation  history,  a 
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mapping  x  on  H  the  six-dimensional  linear  vector  space  spanned  by  components 
e"^  of  the  plastic  strain  tensor  is  introduced 

X  :  H— *  P  (1-4) 

where  P  is  the  positive  class  of  real  numbers.  Other  measures  using  bilinear 
or  nonlinear  maps  involving  both  the  stress  and  plastic  strain  components  are 
in  use.  An  example  is* 

t 

X(<rkl,  e?kl>  =  f  ffkl<T)  de,’kl<T)  <I“B) 

T=-  oo*' 

In  all  cases  the  objective  is  to  define  a  positive  number  k  such  that  it  equals  the 
image  f  under  g  of  boundary  B  of  D.  For  elastic-perfectly  plastic  solids,  k  is 
constant  but,  in  general,  for  stable  material,  k  is  a  monotone  increasing  function 
of  history  of  deformation  and  stress.  In  certain  cases,  the  mapping  g  may  also 
vary  with  plastic  deformation.  This  happens  when  kinematic  constraints  are  pre¬ 
sent.  Theory  of  kinematic  hardening  is  an  instance  in  which  the  reference  point 
in  D,  having  image  zero  in  I,  depends  upon  utress  and  deformation  history. 

Considering,  for  the  present  discussion**  X(e"ij)  =  k, 

g(<rkl)  <  f  =  k  =  *(c"ki)  (1-6) 


*Here  and  in  subsequent  work,  standard  indicial  notation  is  used.  Sum¬ 
mation  on  repeated  indices  is  implied  unless  otherwise  indicated. 

t 

**A  more  general  assumption  uses  X=X(  x  ,  e  "y)  where  k  =  /  *ij(T)  d«"ij(T). 

In  that  case  X  =  .  T  =" 00 

3*  iJ 

ij 
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In  the  Interior  of  D,  c"kl  =  0,  g(o-kl)<  f  =  k  and 

*kl  '  Ekl  <'  'mn> 


=  Ekl  <emn) 

(1-7) 

For  differential  changes  in  stress  and  strain  components,  using  a  superposed  dot 

to  indicate  differential  quantities. 

•  • 

Tkl  =  Eklmn  c’mn 

"  Eklmn  'mn 

(1-8) 

assuming  that  Eki  is  sufficiently  smooth  and  its  derivative  Ekimn,a  tensor  of 

fourth  rank, exists. 

On  the  boundary  B  of  D, 

g(<rkl)  =  f  =  k  =  X  (t"kl) 

(1-9) 

For  g,  X  sufficiently  smooth  in  their  arguments 

*  =  X  <«"kl>  =  hkl'"kl 

d-10) 

g  =  g  ("kl*  *  Okl'ld 

d-11) 

In  the  case  of  elastlo- perfectly  plastic  solids,  hkj  = 

0  and  arbitrary  plastic  straining 

• 

can  occur  for  X  =  0  i.e.  X=  k,  a  positive  constant. 

• 

Also  g  =  f  =  k  requires  g  =  0 

leading  to  the  relationship 

• 

^kl^kl  =0 

(1-12) 

Equation  1-12  requires  the  stress  changes  to  be  in  a  plane  tangent  to  the  hyperplane 

defined  by  g(crkj)  =  k.  However,  for  hw  /  0,for  nonvanishing  c  "kj,  X  >  0  and  g  = 

*  •  •  •  •  #  • 
f  >  0.  This  is  termed  loading  and  g  =  f  =  k  =  X,  For  vanishing  e  "kj,  X  -  0  and 
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once  again  equation  1-12  applies.  This  is  the  ease  of  neutral  loading.  In  all 
eases  g  <0  implies  decreasing  load.  This  is  the  case  of  unloading  and  equation 
1-8  applies  with  =  0. 

1,2.2  Evaluation  of  Incremental  Plastic  Strain  In  Loading 

Equations  I- 10  and  I- 11  suggest  a  relationship  of  the  type 

*"kl  ”  sklmn  'mn 

where  may  depend  upon  c  "mn»  <rmn .  Resolving  emn  into  components  along 

the  boundary  B  and  normal  to  it,  «»e  plastic  strain  is  due  only  to  the  normal  com¬ 
ponent.  Prager  (1949/  showed  that 


independent  of  <rmn  .  Hence  direction  of  c  "jj  Is  independent  of  the  direction  of 
stress  change  given  by  <rmn.  Other  relations  for  plastic  strain  Increments  have 


been  proposed.  Using  a  thermo-dynamic  postulate,  Drucker  (1951)  obtained  the 


normality  rule 


at  g  =  f 


(1-18) 


where  \  is  a  positive  scalar  which  for  rate  Independence  must  be  homogeneous 
of  order  one  in 
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Drucker  used  the  normality  rule  to  evaluate  A .  Defining 


r”  -  T1  r"  r"  V 
e  e  L2  li  ‘  IJ  J 


X  - 


at  g=  f 


(1-16) 

(1-17) 


1  _l£_" 

2  ^  --v 


3<r 


r  11 
2 _ 

IJ. 


Now  defining 


where 


<re  =  ■* 


=  [i  sii  sii]" 


s  =  o-..  .  8,.  ^kk 

^4  1  ''ii  iim  ■  ■  ~ 

3 


Jij  “  uij  "  15 


(1-18) 

(1-19) 


and  writing  _£L  ,  the  slope  of  the  <re,  e  "e  curve  as  H,  equation  1-17  gives,  with 
•  ■■  ^ 


the  normality  rule 


'".i  ■ 


3g'li 


H  * 


■ 


-is. 

3<r 


mn 


.1.6—  1  ^ 


mn 


1  g 

^ij  ^1  Skl 


2  <re  H 


"l 

2 


JJL 


3<ri 


mn 


-UL-l 


i 


(1-20) 


This  formulation  was  used  in  the  so-called  tangent  modulus  methods  e.g.  Swedlow 
and  Yang  (1965). 

Hill  used  the  normality  rule  assuming  X  to  be  a  fourth  rank  tensor  linear 

in  (T^j  and  introduced  a  plastic  potential.  Using  normality  as  well  as  the  condition 
•  •  •  • 

g  =  f  =  k  =  X  on  the  boundary  B  of  D,  Prager  (1949)  obtained 


X  = 


JL 


ILL 


rij 


(1-21) 


'fct" 
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This  formulation  breaks  down  for  elastic-perfectly  solids  where  X  Is  independent 
of  c  "mn.  Felippa  (1966)  obtained  \  in  terms  of  c  ,  increment  in  the  total  strain 
tensor.  In  this  procedure 


where 


Jklmn 


[a*_ 

L*"ii 


-2t.  F 

Bo^ij  'Bo’ij  ‘JM  3<rpq 


r 


• 

emn 

(1-22) 

^4s 

Erskl 

(1-23) 

This  approach  is  valid  for  all  cases  including  perfect  plasticity  and  was  used  by 
Zienkiewicz,  Valliappan  and  King  (1969)  in  developing  finite  element  procedures. 

Using  rate  of  work  equations,  it  is  possible  to  evaluate  X  in  terms  of  stress 
rates  for  materials  of  von  Mises  or  Mohr-Coulomb  type.  Yamada  (1968)  used  this 
approach  for  finite  element  analysis  of  von  Mises  solids.  Using  Drucker  and  Prager's 
generalization  of  Mohr-Coulomb  law,  Reyes  (1965)  developed  the  stress-strain  equa¬ 
tion  for  generalized  Mohr-Coulomb  elastic-perfectly  plastic  solids  under  plane 
strain  conditions.  The  finite  element  procedures  presented  by  Reyes  and  Deere 
(1966),  Baker,  Sandhu  and  Shieh  (1969)  and  those  included  in  Chapter  III  of  this 
report  were  based  on  these  equations.  For  plane  strain 


(1-24) 
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H 

•  b 

’  D11 

D12 

D13  1 

l*11) 

J  °22 

■  ■  2  G 

D21 

D22 

D23 

|  e22  / 

K2 

■  °31 

D32 

D33  - 

(  V12  ) 

where 


Du  =  1  -  1»2  -  2  hi  (Tii  -  h3  <rn2 

E>22  =  1  -  t»2  -  2  hj  o"22  -  I13  «T222 
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4  -  h„  <r 
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and 


D12 

D13 
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h,  = 


°21  =  "  h2  ‘  hl  °22  ‘  hl  "'ll 

D31  =  *hl  " 12  "  h3  " 12  <rll 

D  =  -  h  <r  -h«r  «r 
32  1  12  3  12  22 


3  „  K  -  akk 

2  G 


h3  °22  ffll 


6  3  2^ 


J }  (!  +  9  «2  —  ) 
1  G 


(1-25) 


hn  = 


Tkk 
J2^ 


K  akk 

a  - - r  3q  —  -  - 7“ 

6J,i  G  3J2i 


(1  +  9  a2  J£ ) 
G  ' 


3p  f  K 


J2 


E(l  +  9a2  — ) 
G 


ho  = 


2  J2!  (1  +  9  a2  iL) 
G 


J2  *  i  Slj  S1J 

E,  K,  G  =  elastic  Young's  modulus,  bulk  modulus  and  shear  modulus,  respectively. 


tan  <t> 
a  —  - 


V9+ 12  tan2  $ 


<t>  =  the  angle  of  internal  friction 


1.2.3.  Kinematic  Constraints 

Plane  strain  conditions  impose  a  kinematic  constraint  upon  the  deforming 
solid.  In  relation  to  elastic-plastic  behavior,  a  consequence  is  that  the  yield  sur- 
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surface  from  the  elastic  side  ^nd  plastic  side  do  not,  in  general,  coincide  (Baker 
et.al.1989).  Consider  the  deformation  of  a  body  undergoing  deformation.  F,  the 
set  of  all  admissible  deformation,  is  contained  in  the  six-dimensional  vector 
space  S  spanned  by  components  ry  of  the  strain  tensor.  A  kinematic  constraint 
can  be  written  as 

C  <ckl)  =  0  (1-26) 

and  the  admissible  deformation  is  restricted  to  the  intersection  of  F  with  the  hyper¬ 
plane  in  S  defined  by  equation  1-26.  If  several  constraints  are  present,  the  admiss¬ 
ible  deformation  is  restricted  to 

n 

n  lci  <«kl)  =  °J  (1-27) 


As  the  multiple  Intersection  reduces  the  dimension  of  the  vector  space  by  n,  it  is 
clear  that  n  cannot  exceed  six. 

Consider  a  single  constraint.  In  differential  form  the  equation  is 


ckl  «kl  =  0  (1-28) 

where  coefficients  Cki  depend  upon  tmn. 

As  elastic-plastic  behavior  is  studied  with  reference  to  loading  paths  in 
the  stress  space  V,  it  is  necessary  that  kinematic  constraints  be  rewritten  as 
constraints  on  stress.  Here,  for  no  plastic  strain,  we  simply  use  the  inverse  of 
equation  1-8  to  write  equation  1-27  as 

Ckl  Cklmn  ffmn  =  0  <1-29) 

or  Gmn  ®mn  = 0  <1-30) 
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where 


(1-31) 


^klmn  ”  [  ^klmn] 

Gmn  -  Ckl  Ckhn„  (1-32) 

For  the  case  of  not  all  of  e"^  vanishing,  two  alternative  procedures  are  available. 
Using  definition  of  elastic  strain  tensor, 


Ekl  =  E'kl  +  E  "kl 

If  Tjcjmn  is  the  inverse  of  Sjcjmn  in  equation  1-13 


ckl  *  cklmn  <rran +  Tklmn  fl’mn 

"  £*klpq  +  Tklmn  ®mnpqj  £  pq 

=  K  r  1 

alpq  pq 


(1-33) 

(1-34) 

(1-35) 

(1-36) 


where  1^^  is  a  fourth  rank  identity  tensor.  Thus  the  constraint  is  expressed  by 

C. .  c. .  =  C. .  K. .  c '  (1-37) 

kl  kl  kl  Tclpq  pq 

=  ckl  Kklpq  cpqmn<rmn 


~  Lklmn  rmn  “  0 


(1-38) 
(1-39) 

where  Lklmn  “  ^kl  £*klpq +  ^klrs  ^rspq  j  ^pqmn  (1-40) 

An  alternative  procedure  is  to  use  the  normality  rule  and  to  satisfy  the  constraint 
both  upon  loading  and  unloading  i.  e.  Cy  t  'y  =  0  =  Cy  c  "jj  .  Then  the  first  equa¬ 
tion  is  identical  with  equation  1-29  but  the  second  equation  gives 

(1-41) 

Equation  1-41  may  or  may  not  coincide  with  equation  1-39.  Equations  1-29  and  1-39 
have  linear  relationship  between  incremental  stresses  and  describe  hyperplanes 


Ckl  X  =  ° 

®°kl 
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tangent  to  any  loading  path  in  the  stress  space  V.  As  the  two  equations  are  in 
general  different,  there  is  a  slope  ciscontinuity  in  the  stress  path  as  plastic 
straining  begins  upon  reaching  the  boundary  B  of  D.  We  note  in  particular  that 
proportional  stress  paths  in  V  may  not  be  possible  in  the  presence  of  kinematic 
constraints.  In  the  case  of  linear  elasticity,  let  equation  1-30  define  a  plane 
passing  through  the  origin  in  V.  A  proportional  loading  path  lying  in  this  plane 
is  possible  upto  the  point  of  intersection  with  boundary  B.  Beyond  that,  upon 
loading,  stress  path  has  to  be  in  the  surface  determined  by  Equation  1-39  and 
this  will  in  general  be  non-planar.  If  loading  is  continued  to  a  certain  point 
along  this  surface,  unloading  therefore  will  be  along  a  path  lying  in  a  plane 
parallel  to  the  original  loading  plane  but  different  from  it  and  not  passing  through 
the  origin.  Thus  unloading  to  initial  state  is  impossible.  This  corresponds  to 
setting  up  of  residual  stresses  corresponding  to  kinematic  constraints. 

Specifically  considering  plane  strain  conditions  and  elastlc-perfectly 

plastic  Mohr-Coulomb  material,  the  mapping  g  from  the  set  of  all  admissible 

stress  states  to  the  positive  Interval  £o,  f  ]is 

«Tld{  i 

g  (o-ij)  =  Q  3  +  J22  (1-42) 

«Tkk2 

J2  =  i  ‘’’ij  <rij  “  ”g”  <I_43> 

Linear  isotropic  elasticity  implies,  for  ?"|j  =  0 

eij  =  c'ij  =  ~Fg  <rij  '  “TT^  *kk  (I"44) 

Plane  strain  condition  implies 

c13  ”  c23  =  c33  ~  0  (I_45) 
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From  equations  1-44  and  1-45 


I 
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I 
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I 


(1-46) 

(1-47) 


(1-48) 


Equations  1-47  and  1-48  define  different  surfaces  in  a  three-dimensional  space 
spanned  by  or11#  e^,  «,33.  Let  their  intersections  with  B  be  respectivelyy  P  and 
Q.  The  stress  path  is  constrained  to  lie  in  the  plane  defined  by  equation  1-47  for 
stress  states  in  the  interior  of  D  and  for  neutral  loading.  For  plastic  deformation 
to  occur  the  stress  path  must  lie  in  the  surface  defined  by  equation  1-48.  For  a 
continuous  stress  path  to  be  possible,  P  and  Q  must  coincide.  In  general  this 
is  not  the  case.  Figure  1-7  illustrates  the  difference  between  the  surfaces  P 
and  Q  for  Mohr-Coulomb  plane  strain  case. 

Considering  that  the  stresses  «r33  do  not  contribute  to  energy /work  of  the 
system,  it  appears  reasonable  to  assume  that  progress  from  P  to  Q  is  possible 
with  gradually  increasing  the  value  of  <r33.  This  ^ould  amount  to  following  the 
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from  Plastic  Side 


boundary  B.  Growth  of  0-33  in  progress  from  P  to  Q  and  the  behavior  upon  unloading 
are  not  clearly  understood.  Later  investigations  may  throw  light  on  this  aspect  of 
material  behavior.  For  the  purpose  of  the  computer  program  in  Chapter  in, it  is 
assumed  that  elastic  loading  in  plane  strain  can  be  continued  upto  a  point  from  which 
the  transition  to  plastic  plane  strain  loading  is  possible  merely  by  adding  a  residual 
value  of  <r33.  Referring  to  Drucker  and  Prager  (1952),  this  is  given  by 


k 


f  g  (<rll»  <r22,T12^ 

[2 
v  11  ~  <r22 
. -"T ~ 


(1-49) 


1. 3.  Stress-Strain  Behavior  of  Jointed/Cracked  Rock 

Mathematical  simulation  of  behavior  of  jointed  rock  must  allow  for  closing 
of  pre-existing  open  joints  under  compressive  loads  followed  by  linear  elastic 
behavior  upto  initiation  of  fracture.  After  fracture  occurs,  the  material  cannot 
take  any  tension  locally  in  the  direction  normal  to  the  plane  of  crack.  Non-monotonic 
loads  may  involve  closing-cracking-closing  cycles. 

The  finite  element  method  has  been  applied  to  jointed  rock  (Anderson  and 
Dodd,  1966;  Goodman  et  al,  1968;  Duncan  and  Goodman,  1968;  Malina,  1971). 
Anderson  and  Dodd  used  pin  ended  one-dimensional  elements  across  a  fault  to 
allow  compressive  stresses  to  be  transferred  in  the  direction  normal  to  the  fault. 

The  fault  plane  was  assumed  to  have  no  resistance  against  shear  or  tensile  loads. 
This  capability  is  now  routinely  incorporated  in  most  finite  element  programs.  A 
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two-dimensional  'soft'  material  element  has  long  been  used  to  represent  weak 
joint  planes  in  rock.  Duncan  and  Goodman  (1968)  object  to  this  on  the  basis  of 
large  number  of  elements  needed  to  ensure  a  reasonable  'aspect  ratio'  in  shnpe 
of  elements.  This  becomes  a  problem  for  elements  representing  very  thin  joints. 
A  one-dimensional  element  with  shear  and  normal  stiffness  characteristics  was 
proposed  by  Goodman  et  al  (1968)  to  eliminate  this  objection.  Recently  (1971), 
the  same  investigators  have  introduced  nonlinear  properties  in  this  type  of  ele¬ 
ment.  This  approach  is  quite  effective  for  the  case  of  pre-existing  joints  in  rock. 
For  well  defined  orthogonal  joint  systems,  an  orthotropic  continuum  approach 
was  suggested  by  Duncan  and  Goodman  (1968).  Christian  is  credited  (Einstein, 
Bruhn  and  Hirschfeld,  1970)  with  development  of  an  element  capable  of  simulating 
constant  shear  and  residual  shear  characteristics. 

In  all  these  investigations,  a  distinct  set  of  elements  is  used  to  represent 
the  joint.  This  is  alright  for  pre-existing  joints  but  is  impracticable  for  dis¬ 
continuities  arising  as  a  result  of  fracture  under  applied  load.  To  use  the  same 
procedure  both  for  pre-existing  joints  as  well  as  post  failure  cracks,  it  is  nec¬ 
essary  to  allow  cracks  and  joints  within  elements.  Then,  the  mesh  layout  is 
more  flexible  and  arbitrary  failure  laws  can  be  used.  Malina  (1971)  used  this 
approach  to  study  failure  along  joint  planes  and  then  went  on  to  compute  the 
amount  of  slip  and  accompanying  stress  redistribution  on  the  basis  of  deforma¬ 
tion  or  slip  theory  of  plasticity. 
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Apparently,  a  bimodular  analysis  procedure  (Sandhu  and  Wilson,  1970) 
can  be  used  to  represent  pre-existing  joints  as  well  as  fractures.  Bimodularity 
would  be  dependent  upon  the  joint  opening.  However,  noting  that  fractured  or 
open  jointed  rock  has  no  resistance  to  tension  in  the  direction  normal  to  that  of 
fracture,  a  simple  approach  following  the  procedure  introduced  by  Zienkiewicz 
et  al  (1968)  is  more  convenient.  The  'no  tension'  method  of  Zienklewicz  consists 
of  first  obtaining  a  solution  assuming  the  system  to  be  linear  elastic.  Then  the 
elements  in  tension  are  relieved  of  the  tensile  stresses  by  application  of  self- 
equilibrating  forces  in  elements  and  at  nodal  points.  This  gives  an  Iterative 
scheme  for  redistribution  of  loads  to  surrounding  rock  and  a  lower  bound  to 
the  exact  solution.  This  approach  is  essentially  an  orthotropic  continuum 
approach  with  the  orthotropy  being  applied  to  individual  elements  depending 
upon  the  orientation  of  the  fracture  plane.  The  fracture  plane  defines  also  a 
plane  of  material  orthotropy.  The  relationship  between  principal  stresses  and 
strains  can  be  written  as 
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(1-50) 


or  symbolically 

<rp  =  Cp  Cp  (1-51) 

The  laws  of  transformation  of  stress  and  strain  give, 
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and 
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where  v%t  <Ty,  are  components  of  stress  and  e^,  ty,  y^  are  oomponenta  of 
strain  In  x,  y  coordinate  system  and  $  is  the  angle  between  the  prlnoipal  direction 
1  and  x-axls.  Symbolically,  the  above  equations  are 


«r  =  JT  e- 


*P-  J  ' 


(1-54) 

(1-55) 


Substitution  in  1-61  gives 


v  -  jT  cp  J  c 


(1-55) 


Cc 


where 


C  *  J1  C_  J 


(1-57) 

p  -  (1-58) 

Equation  1-58  gives  the  transformation  for  stress-strain  relation  for  principal 
direction  to  any  arbitrary  choice  of  coordinates.  The  matrix  C  la  singular  only 
for  0  =  0  or  90*.  It  is  thus  possible  to  use  the  relationship  In  principal  stresses 
and  principal  strains  as  the  starting  point. 

In  finite  element  analysis  procedures,  the  stiffness  mxtrlx  for  the  system 
Is  the  sum  of  element  stiffnesses. 


M 

K  =  Z  km 
m  =  l 


(1-50) 


23 


where  K  is  the  system  stiffness,  km  is  the  stiffness  of  the  nth  element,  and  £ 
is  viewed  as  a  direct  stiffness  summation  operator.  Further,  element  stiffness 
is  related  to  constitutive  relationship  through  the  equation 


„m 


,/[ 

'm 


bT  c  b]  dV 


(1-60) 


where  b  is  the  matrix  relating  strains  to  nodal  point  displacements,  and  V  is  the 
volume  of  the  element.  Using  Equation  1-58,  the  integrand  in  1-60  can  be  written 


bTC  b  =  bT  JT  Cp  J  b  (1-61) 

=  bT  JT  o-p  (1-62) 

=  BT  <rp  (1-63) 

T  T  T 

where  B  =  b  J  relates  principal  stresses  to  nodal  point  forces. 

Occurrence  of  fracture  in  an  element  reduces  its  ability  to  take  tensile 
stresses  normal  to  the  fracture  plane.  Also  there  can  be  no  shear  transmitted 
across  a  crack  which  therefore  is  a  principal  direction.  Thus,  it  is  reasonable 
to  reanalyze  the  system  assigning  an  orthotropic  constitutive  relationship  and  a 
presecrlbed  principal  direction  to  the  element  containing  a  fracture.  The  proce¬ 
dure  is  to  be  repeated  until  no  further  fracturing  occurs  under  a  given  load.  To 
allow  for  nonlinearity  introduced  by  progressive  cracking,  incremental  procedures 
are  required. 

In  Zienkiewicz  et  al  (1968),  a  change  in  element  stiffness  was  considered 
equivalent  to  a  pseudo-load.  Thus  an  iterative  solution  scheme  was  set  up  in 
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which  each  iteration  only  involved  a  back- substitution  operation.  The  pseudo¬ 
loads  were  computed  as  equivalent  to  tensile  principal  stresses.  This  is  satis¬ 
factory  when  both  principal  stresses  are  tensile.  However,  when  only  one  of 
the  principal  stresses  is  tensile,  use  of  pseudo-load  corresponding  to  one  prin¬ 
cipal  stress  introduces  a  non- symmetric  constitutive  law.  Actually  if  the  phy¬ 
sical  concept  of  'unloading  without  any  displacements'  be  followed,  a  change  In 
the  second  principal  stress  corresponding  to  Poission's  effect  due  to  the  first 
stress  must  be  Included.  This  modification  is  included  in  the  computer  program 
in  Chapter  IV. 
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CHAPTER  n 

THE  FINITE  ELEMENT  METHOD 


Chapter  II.  The  Finite  Element  Method 
2.1.  Basic  Concepts 

A  boundary  value  problem  can  be  stated  in  the  form 

A  u  =  f  on  F  (II-l) 

where  u  is  the  unknown  function  to  be  determined,  A  is  an  operator,  and  f  is  the 
'forcing'  function.  F  is  the  domain  of  interest  and  may  be  an  open,  connected, 
bounded  spatial  region  embedded  in  R3  or  in  a  cartesian  product,  R3  x  £o ,  <») 
where  £o,  <»)  is  the  non-negative  time  interval.  In  addition  to  the  field  equation 
H-l,  there  will  be  some  conditions  to  be  satsifed  on  boundary  S  of  F.  For  A 
linear  positive,  it  can  be  shown  that  equation  n-1  has  a  unique  solution.  Nec¬ 
essarily,  any  approximate  solution  will  in  general  not  coincide  with  the  unique 
solution  of  II-l  and  consequently  no  approximate  solution  is  expected  to  satisfy 
the  field  equation  as  well  as  the  boundary  conditions  completely. 

Solutions  to  engineering  problems  as  well  as  the  forcing  functions  are 
in  general  bounded  and  therefore  belong  to  L2  ,the  space  of  square  lntegrable 
function.  However  u  may  be  contained  in  a  subset  D  of  L2  such  that  A  is  defined 
on  D.  We  assume  that  D  is  dense  in  L2.  If  the  set  of  functions  |  k~  1*2f  «j 
is  an  orthonormal  basis  in  D,  then  any  function  u  can  be  expressed  as  an  infinite 
sum;  N 

u  =  2  ak  <fe  01-2) 

k  =  1 

A  scheme  to  generate  approximate  solutions  is  to  use  a  finite  set  of  terms  in  the 
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<n-3) 


infinite  sum  above.  Thus,  we  use 

N 


as  an  approximation.  The  approximation  process  then  consists  of  appropriate 
choice  of  N,  ^  and  the  coefficients  ak,  Several  alternative  procedures  are  avail¬ 
able.  The  finite  element  method  is  a  special  process  of  selection  of  a  finite  sub¬ 
set  of  the  basis  |  .  The  coefficients  a^  are  evaluated  by  an  extension  of  Ritz 

method  or  other  standard  procedures. 

The  finite  element  method  is  well  documented  in  literature  (Zienkiewicz, 
1967;  Bell  and  Holand,  WPAFB Conference,  1965,1968;  Felippa,  1966;  Clough,  1960, 
1965).  Its  theoretical  basis  (Oden,  1969;  de  Arantes  e  Oliveira,  1968,  Zlamal, 
1968;  Melkes,  1970)  and  relationship  to  variational  principles  (Melosh,  1963,  Plan 
and  Tong,  1969)  have  been  examined.  Essentially,  a  finite  element  idealization 
partitions  the  spatial  region  F  into  a  finite  number  of  nontrivial  discrete  elements 
or  subregions.  The  geometry  of  the  elements  is  defined  by  a  set  of  points  in  space 
called  the  nodal  points  of  the  system. 

Over  an  element  e  let  an  approximation  to  u  be 

Ne 

ue  =  Z  ake  0ke  (n-3) 

k  =  1 


or  in  matrix  form 


©  T  — 

where  { <t>  }  is  a  row  vector  consisting  of  as  its  elements  and  {a6} is  a 


“e  *  [♦*  f  {  a6} 


(H-4) 
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column  vector  of  coefficients  a^®.  Evaluating  the  function  at  nodal  points 

)  “t®}  =  j  «tc  j 1  /  a®  }  <n-S) 

where  |  Uj®  |  is  the  vector  of  nodal  point  values  of  the  function  and|"$ie  J  ^  is  the 

0  T 

matrix  of  base  functions  evaluated  at  each  nodal  point.  Rows  and  columns  °f[$i  ] 
are  linearly  independent.  If  square,  the  matrix  is  Invertible,  If  the  number  of 
nodal  points  is  not  equal  to  the  number  of  independent  base  functions,  a  least 
square  procedure  can  be  used  for  inversion.  Hence,  we  can  write 

[my  M 

-  [A]-1  {  Uie|  <n-6) 

where  A  =  £<^eJ  T 

Substituting  n-6  in  11*5 

T 

»®  =  { ?e }  [a]  _1  |  U,®  I  <n-7) 

=  {*T{Ul*}  <n"8) 

where |  (/)°  |  can  now  be  regarded  as  a  set  of  Interpolating  functions  relating  nodal 
point  values  of  a  function  to  the  value  of  an  arbitrary  point  within  the  element. 

2.2.  A  Potential  Energy  Formulation 

Wo  nesurnc  the  rock  continuum  or  'dlscontlnuum'  to  be  stepwise  linear  for 
nufflotontly  mull  stops  In  loading.  For  such  a  case  the  governing  equations  are 

(II-9) 

(II- 10) 

<n-ll) 


ki,k 
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kl 
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*lk  *  p  F1  *  0 
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where  or^ ,  ,  Ej^i  *  Fi  »  ui  are  components  .respectively, of  the  symmetric 

stress  tensor,  the  symmetric  strain  tensor,  the  isothermal  elasticity  tensor,  the 
body  forces  vector  per  unit  mass,  and  the  displacement  vector,  p  is  the  mass 
density  and  5^  is  the  kronecker  delta,  tr^  are  components  of  Initial  stress  cor¬ 
responding  to  zero  displacement  and  ir  is  the  pore  pressure.  Potential  energy 
formulation  uses  the  functional 


»■/[■ 


ijEijklekl  “  ui<rij,J“2tiJorij  +  ui,<rij 

-2ujpFj  -  2  Uj  7rti+  2  ejj  ryj  dF 


+ 


A 

ui  (^ij  nj  -  2  tj  )  d  s 


A 

(Ui  -  2  Ui ) 


aijUj 


ds 


(II- 12) 


where  we  have  included  the  boundary  condition 

OTjj  n^  =  on  8^^  (11-13) 

Uj  =  Uj  on  s2  (II- 14) 


s1(  s2  are  complementary  subsets  of  S  the  boundary  of  F  and  nj  are  components 
of  unit  vector  normal  to  surface. 

Symmetry  of  the  field  equations  leads  to  the  functional 


ft  = 


/[' 


-  2 


Hj  Eijki  *kl  2  eij  ,rij  +  2  'H.j  ^ij  “  2  'H  p  Fi 
-2u1h->1  +  2cij7iJ]dF 

f  A  f  A 

f  Ui  t|  d  s  -  2  I  (Ui  -  Ui  )  <Tij  nj  d  s 
'1  B2T 


(n-15) 


Further  assuming  that  we  restrict  our  choice  of  ey  ,  14  such  that  II- 11  and  11-14 
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arc  identically  satisfied,  the  functional  in  11—15  reduces  to 

Eijkl  %1  -  2uiP  Fi  -  2Ui  7rfif  2  Uijffy  J 


dF 


-  2 


/m. 


ds 


(n-i6) 


Replacing 


f  by  m  =  1  /* 


«fm  where  Fm  represents  the  subregion  or  element 
m,  and  using  suitable  interpolation  scheme  to  express  the  integrand  in  terms  of  nodal 
point  vectors  of  displacement,  vanishing  of  variation  of  the  functional  yields  the  matrix 


equation 
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Components  of  element  stiffness  matrix  and  load  vectors  are: 
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i  0  i  T 

and  are  components  of  a  matrix  formed  by  the  row  vectors  |  <t>  |  corresponding 
to  each  degree  of  displacement  freedom.  The  vectors  |Lm} ,  |Pm| ,  |Qm } ,  |Tm| 
represent  the  contribution  to  the  load  vector  made  respectively  by  the  body  forces, 
the  pore  pressure  gradients,  the  initial  stresses  and  the  boundary  loads  in  the  ele¬ 
ment  m. 


2 . 3  Incremental  Analysis 

In  case  of  Incremental  construction  and  vncremental  application  of  loads, 
the  loads,  stresses  and  displacements  for  any  incremental  step  can  be  written  as 
|Rn|  *  {^n}  and  |un|  •  Then  for  the  next  stage,  |irn|  and  junj  can  be  regarded 
as  the  Initial  stresses  and  the  initial  displacements  for  the  structural  system.  Thus 
the  matrix  equations  are 


where 


(Xu]  {Vl  -  =  {AHn}  <n-25' 

ml  1 

<n-26) 


m 

ARn  "  S 
m  =  l 


ALm  +  APm: 


AQ1"  +  AT™ 


and  |ALm^  ,  |APm|,|AQm|,  j  ATm|  are  increments  in  the  respective  quantities. 

For  elastic-plastic  analysis,  the  stiffness  depends  upon  stress  and  has  to  be 
re-evaluated  at  small  increments  of  load.  To  ensure  manageable  computation,  the 
increments  are  kept  at  the  largest  practicable  without  loss  of  accuracy. 
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CHAPTER  m 


COMPUTER  PROGRAM  FOR  PLANE  STRAIN  ANALYSIS  OF 
ELASTIC -PLASTIC  MOHR-COULOMB  MATERIALS 


Chapter  III.  Consput**’  K  -ig*  m  "  Plan*  Strain  Analysis  of 
Elastic- Plastic  Molir-^oulomb  Materials 

3.1.  Organization 

Computer  program  described  hero  is  based  on  the  theory  presented  earlier 
in  this  report.  The  program  is  written  in  Fortran  IV  language. 

The  program  is  intended  to  calculate  stresses  and  strains  for  a  plane  strain 
problem  in  rock  mechanics.  Mohr-Coulomb  yield  criterion  has  been  used.  It 
makes  allowance  for  the  boundary  conditions,  residual  stresses,  stresses  due  to 
temperature  change,  and  varying  pressure  boundaries.  The  structure  may  consist 
of  different  materials.  It  uses  Wilson's  (1965)  quadrilateral  elements  and  generates 
stiffness  in  line  with  the  integration  procedures  discussed  by  Felippa  (1966). 

The  principal  program  called  MAIN  controls  all  the  data  input  and  control 
information.  It  does  the  basic  system  initialization  and  prints  the  control  data  and 
material  and  geometrical  properties  of  the  structure.  Stiffness  formulation,  equa¬ 
tion  solving  and  stress  calculations  are  done  by  the  subroutines  called  by  MAIN. 
3.12.  Stiffness  and  Load  Matrices 

Stiffness  matrix  for  each  analysis  is  computed  in  blocks  by  the  subroutine 
STIFF.  For  the  element  stiffness  it  calls  QUAD  for  triangular  and  quadrilateral 
elements  which  have  been  allowed  by  this  program.  The  element  stiffness  is  added 
to  the  total  stiffness  using  the  direct  stiffness  technique.  Concentrated  loads  are 
included  in  the  load  matrix.  Equations  are  modified  for  the  displacement  boundary 
conditions  by  calling  subroutine  MODIFY.  For  the  stress-strain  matrix  QUAD 
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calls  8TUSTR.  With  the  constitutive  iaw  being  available,  stiffness  of  the  system 
ami  load  matrix  arc  computed  in  die  S'l  iFi  subroutine. 

3. 13.  Calculations  of  Displacements 

After  the  stillness  and  load  matrices  for  a  stage  have  been  computed,  the 
resulting  equations  arc  solved  by  c  tiling  subroutine  BANSOL.  This  uses  Gaussian 
elimination  technique  for  banded  equations  by  Wilson  (1963).  In  this  the  trlangulari- 
?  at  ion  of  stiffness  matrix  is  done.  Back  substitution  through  the  triangulized  matrix 
gives  the  solution. 

14.  Calculations  <>r  Stresses 

In  the  <irst  cycle  i  purely  >‘138110  solution  is  obtained  for  the  problem  by  solving 

[  K°]  { r  (  {K}  <M-1) 

who?  *;[k'  Jis  the  olnr.ttc  stiffness 

j  r  }  is  the  ills  placement  vector 
JR  }  is  the  load  vector. 

This  can  be  done  easily  by  assuming  all  the  elements  to  be  elastic  to  begin  with. 

As  the  problem  is  a  nonlinear  one,  therefore  this  solution  will  not  be  correct.  In 
our  analysis  the  system  is  assumed  to  be  stepwise  linear  between  the  yielding  of 
one  element  to  the  other.  This  is  assumed  not  to  cause  any  significant  error.  In 
Fig.  3-1  point  A  represents  the  initial  stresses  and  C  the  final  stresses  in  an  cle¬ 
ment.  The  curve  f  k  represents  the  yield  surface.  For  those  elements  which 
become  plastic  under  this  loading  AC  meets  the  nurface  f  -  k  at  pt.  B.  It  is  seen 
easily  that  for  the  element  It  is  not  possible  to  be  loaded  to  point  C  but  it  can  Ik* 
loaded  to  point  B  only,  assuming  proportional  loading. 
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I.ct  S,  ,  ■ .  .  1h  called  the  HtreBH  ratio.  To  calculate  S„  we 

r  |AC|  r  ' 

know  that 

AC  =  (o-ij),  -  (<r1j)i 
OB  ^  oli  ‘  A  B 

=  OA  4  Sr  •  A<* 

+  <«rij)1  *  sr  <°*ij  >f  - 

*ij  =  <Vj  <l"Sr>  *  Sr  >  f 
As  the  point  B  lies  on  the  yield  surface  f  -  k 

f  [<Vi<l’Sr>  "  Wf]"  k  (ID‘2) 

From  equation  (III— 2)  the  value  of  Sr  can  be  calculated.  This  stress  ratio  repre¬ 
sents  the  fraction  by  which  the  increment  in  stress  is  to  be  scaled  to  bring  the  final 
load  on  the  yield  surface. 

Value  of  Sr  is  calculated  for  all  the  elements,  The  element  in  which  the 
final  stress  state  is  farthest  from  the  yield  surface  will  have  the  minimum  stress 
ratio.  If  we  scale  down  the  displacements  and  stresses  in  this  ratio,  we  shall  have 
the  stresses  and  strains  precisely  at  the  point  when  the  system  has  its  first  element 
just  going  into  the  plastic  region  from  purely  clastic  system. 

In  the  next  step  the  element  having  stress  ratio  equal  to  the  minimum  value 
is  assumed  to  be  plastic.  To  economize  on  computer  time  all  such  elements  which 
have  their  value  of  stress  ratio  in  the  vicinity  of  minimum  were  also  allowed  to  go 
plastic.  As  the  stress-strain  matrix  is  known  the  stiffness  is  calculated  again  and 

NH=W 


equation 


(in-:*) 


is  solved.  This  procedure  is  repeated  until  the  whole  load  has  been  applied  to 
the  system  and  cumulative  stresses  and  displacements  calculated.  The  stresses 
in  (i-l)th  step  become  the  initial  stresses  for  ith  step. 
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3.2.  Input  Data  Preparation 


1.  Control  Card  (A6).  This  card  will  carry  the  characters  START  in  columns 
1-5.  This  will  start  the  processing  of  the  data  deck  which  consists  of  the 
following  set  of  cards. 

2.  Job  Title  (72 H).  This  card  will  give  the  descriptive  identification  for 
the  job. 

3.  Control  Information  (415,  2F10.2,  15) 


Information  Columns 

Total  number  of  nodal  points  1-5 

Total  number  of  elements  6-10 

Number  of  different  materials  11-15 

Number  of  pressure  boundary  cards  16  -  20 

Body  Force  in  x-direction  21-30 

Body  Force  in  y-direction  31  -  40 

Number  of  Approximations  41  -  45 


4.  Material  Property  Cards.  One  set  of  2  cards  is  provided  for  each  material. 
In  each  set: 

a.  first  card  (115,  F10.0)  will  give  the  following  information 

Material  identification  number  1-5 

Maas  density  of  the  material  6-15 

b.  The  second  card  will  carry  the  following  information  (4F10.0) 


I 
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Information 


Columns 


Elastic  Modulus  1-10 

Poisson's  Ratio  11  -  20 

Cohesion  21-30 

Friction  Angle  in  Degrees  31-40 


5.  Nodal  Point  Cards  (15,  F5.0,  5F  10.0). 

One  card  for  each  nodal  point  with  the  following  information! 


Nodal  Point  number 

1-5 

Type  of  Nodal  point 

6-10 

X-ordinate 

11  -  20 

Y-ordinate 

21  -  30 

XR 

31  -  40 

XZ 

41  -  50 

If  the  number  in  columns  6  -  10  is 

Zero  XR  is  the  specified  X-load  and  XZ  is  the  specified  Y-load 

1  XR  is  the  specified  X-displacement  and  XZ  is  the  specified  Y-load 

2  XR  is  the  specified  X-load  and  XZ  is  the  specified  Y-displacement 

3  XR  is  the  specified  X-displacement  and  XZ  is  the  specified  Y-displacemc 

All  loads  are  considered  to  he  total  forces  acting  on  an  element  of  unit  thickness. 

Nodal  point  cards  must  be  in  numerical  sequence.  If  cards  are  omitted,  the 

omitted  nodal  points  are  generated  at  equal  intervals  along  a  straight  line  between 
the  defined  nodal  points.  The  type  of  the  nodal  point,  as  well  as  XR,  XZ,  are  set 
equal  to  zero. 
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6. 


Element  Cards  (615).  One  card  for  each  element  will  provide  the  following  data. 


Information  Columns 

Number  of  element  1-5 

Nodal  point  I  6-10 

Nodal  point  J  11-15 

Nodal  point  K  16-20 

Nodal  point  L  21-25 

Material  type  26  -  30 


Nodal  points  I,  J,K,L  arc  comers  of  each  individual  element  in  a  counter¬ 
clockwise  order  for  a  right  handed  system  of  coordinates.  For  triangular 
elements  set  nodal  point  L  same  as  nodal  point  K.  The  element  cards  must 
be  in  the  numerical  sequence.  Any  cards  that  are  omitted  will  be  automatically 
generated  in  the  program  by  incrementing  each  of  the  It  JtKf  and  L  nodal  points 
by  one.  The  material  type  wil!  be  taken  the  same  as  for  the  last  element  defined. 

7 .  Pressure  Boundary  Cards  (2I5y  2F10 . 0) .  One  card  for  each  boundary  element 
which  is  subjected  to  a  normal  pressure  will  carry  the  following  information: 


Information  Columns 

Nodal  Point  I  1-5 

Nodal  Point  J  6-10 

Normal  Pressure  at  I  11  -  20 

Normal  Pressure  at  J  21-30 
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As  shown  in  the  sketch,  the  boundary  element  must  be  on  the  left  as  one 
progresses  from  I  to  J.  Surface  tensile  forces  is  input  as  a  negative  pressure. 

Output  Information: 

The  following  information  is  developed  and  printed  by  the  program: 

1 .  Reprint  of  input  data 

2.  Nodal  point  displacements 

3.  Stresses  at  the  center  of  each  clement 
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o  o 


3.3  Program  listing 

PLANE  STRAIN  ANALYSE  Of  ELA  STIC  -  PLASTIC  MOHR-COULOMB  MATERIALS 
IMPLICIT  REAL*8( A-H.O-Z ) 

COMMON  ACELR, ACELZ , VOL , TEMP,  SIGI C 400, 71 VHE0(18) ,E(4, 12) *R0I 121  t 
♦R(500),Z(500) »URI 500), UZ I  500) ,CODE( 5001 ,T< 5001. PRI 100.21 » 

♦NUMMAT  ,  NUMPC * N, MT YP E ,KKK, NUMNP, NUMEL  »NNN* 

*  I BC ( 100 1  *  JBC  (  1 00 ) » MT AG( 400  I 

COMMON/ ARG/  RR I  5 )  ,  ZZ  1 5  I ,  SUO ,  10 )  ,  PI  8  ) ,  ST  I  3, 10)  ,C  I  3 , 3  I ,  SIG17)  ,  EE  1 41 
♦,BOtlOOO»,SRl,SR2, 

♦  RAT  1 0(400 ) ,LM(4) , I X 1 400 » 5  I , XC , YC 
COMMON  /9ANARG/  MBANO, NUM8LK » B(  108  I , A<  108 , 54 1 
DIMENSION  WORD( 2 ) 

DATA  WORD/  6HSTART  ,6HST0P  / 

CALL  ERRSET( 20 7,256,-1*1) 

CALL  ERRSET ( 208 , 256, -1 , 1 ) 

C 

C 

5  READ  (5*1006)  W0R01 

IF  ( WORDl„EQ.WORO( 1 1  I  GO  TO  50 
IF  ( W0R01.EQ.WORD(2) I  STOP 
GO  TO  5 

50  READ  (5*1000)  HED, NUMNP, NUMEL , NUMMAT .NUMPC »ACELR* ACELZ *NP 
WRITE (6,2000)  HED, NUMNP, NUMEL, NUMMAT .NUMPC ,ACELR, ACELZ ,NP 
DO  55  M*l, NUMMAT 
READ  (5,1001)  MTVPE.ROl MTYPE I 
WR I T E (6 • 200 1 )  MTYPE, RO(MTYPE) 

READ  (5,1002)  ( E ( J , MTYPE ) , J* 1 , 4) 

WR I TE (6,2002 )  ( E ( J .MTYPE ) , J«1 ,4) 

55  CONTINUE 
C 

WRITE  (6,2003) 

L*0 

C 

60  READ  (5,1003)  N,COOE(N) ,R(N),Z(N),UR(N),UZ(N),T(N) 

NL*LM 

ZX*N-L 

DR* ( R ( N l-R ( L l )/ZX 
DZ=(Z(N)-Z(L) )/ZX 
DT*(T(N)-T(LI l/ZX 
70  L*L*l 

IF(N-L)  100,90,80 
80  CODE ( L ) *0,0 
R ( L ) *R ( L-l I  FOR 
Z(LI *Z(l-l |40Z 
T(L)*T(L-1) +DT 
UR(LI*0.0 
UZ(L)*0.0 
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GO  TO  70 

90  WRITE  (6,2004)  ( K ,COOE ( K ) , R ( K I , Z ( K) ,UR< K I ,UZ ( K ) * T( K ) ,K-NL , Nl 
I F ( NUMNP-N )  100,110,60 
100  WRITE  (6,2005)  N 
CALL  EXIT 
110  CONTINUE 

WRITE  (6,2006) 

N=0 

130  READ  (5,1004)  M, ( I  XI M, I ) , I  =  l , 5) > ( SI GI ( M, I ) , I =  1 ,4) 

ZX=M-N 

DO  135  1*1,4 

135  SIG( I )*( SIG1(M,1)-SIGI(N,I ) )/ZX 
140  N=N*l 

IF  (M-N)  170,170,150 
150  IX(N,l)=IX(N-l,l)M 
IX(N,2)=IX(N-I,2)*l 
I  X  (  N , 3 ) *  I X( N-l  »3)4,l 
I  X  (  N  , 4  )*  1  X(  N- 1 , 4 )  ♦  l 
IX( N,5)*IX(N-l,5) 

DO  160  1=1,4 

160  SIGI(N,  I  )=SIGI(N-1,I)+SIG(  I  ) 

170  WR I TE ( 6 , 2007 )  N, ( I X( N , I ) , I = l , 5 ) , ( S IGI ( N, I ) , I *1,4) 

IF  (M-N)  180,180,140 
180  IF  (NUMEL-N)  190,190,130 
190  CONTINUE 

IF  (NUMPC)  290,310,290 
290  WRITE  (6,2008) 

DO  300  L* l, NUMPC 

READ  (5,1005)  IBC( L ) , JBCl L ) • PR ( L , 1 ) , PR( L  *  2 ) 

300  WRITE  (6,2009)  I8C( L ) , JBC ( L ) ,PR( L  « 1 ) ,PR( L ,2 ) 

310  CONTINUE 

J=0 

DO  340  N*l,NUMEL 
MTAG(N) *0 
SIGI (N,5)*0. 

SIGI(N,6)*0. 

SIGI (N,7)*0. 

DO  340  1=1,4 
DO  325  L=l,4 
Kk= I ABS (  I  X ( N ,  I  ) - 1  X ( N , L ) » 

IF  (KK-J)  325,325,320 
320  J*KK 
325  CONTINUE 
340  CONTINUE 

MBAND=2*J*2 
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WRITE (6, 1007)  MBANO 
00  350  N* 1 » NUMNP 
B0( 2*N- 1  )*0  • 

350  B0!2*N)=0. 

C 

SR  1*  L  .0 
SR2=0.0 

00  500  NNN*1,NP 
KKK=0 

CALL  STIFF 
CALL  BANSOL 
CALL  STRESS 
00  400  N* l, NUMNP 
NN=2*N 

BO(NN-1)*BO(NN-1  M-B(NN-l) 

BO(NN|=BO( NN) *B!NN) 

400  CONTINUE 

WRITE! 6, 20 10)  (N,B0!2*N-1) ,B0!2*N), N«l, NUMNP) 

IF(KKK.EQ.O)  CALL  EXIT 
500  CONTINUE 
GO  TO  5 
C 

1000  FORMAT  l 18A4/4I5,2F10.2, 21 5) 

1001  FORMAT  ( 1 1 5 *  IF  10 *0 ) 

1002  FORMAT  16F10.0) 

1003  FORMAT  (1 5 , F5.0, 5FI0.0 ) 

1004  FORMAT! 61 5,4F10.0 ) 

1005  FORMAT  (2I5.2F10.0) 

1006  FORMAT ( A6 ) 

1007  FORMAT!  •  BAND  WIDTH  FOR  THIS  DATA  »  •  ,  15  » 

2000  FORMAT  ! IHl  IBA4/ 

1  30H0  NUMBER  OF  NODAL  POINTS - 13  / 

2  30H0  NUMBER  OF  ELEMENTS - 13  / 

3  30H0  NUMBER  OF  DIFF.  MATERIALS—  13  / 

4  30H0  NUMBER  OF  PRESSURE  CARDS—  13  / 

5  30H0  X-ACCEL ERATION - F12.4/ 

6  30H0  Y-ACCELERATION - E12.4/ 

7  30H0  NUMBER  OF  APPROXIMATIONS—  I12» 

2001  FORMAT  ( 17H0MATER IAL  NUMBER  I  I3«  15H,  MASS  DENSITY#  E12.4) 

2002  FORMAT! 16H0ELASTIC  MOOULUS  14X  2HNU  BX  8HC0HESI0N  2X  14HFRICTI0N  A 
♦NGLE/I2E16.5.2F16.5M 

2003  FORMAT  (UIH1N0DAL  POINT  TYPE  X  ORDINATE  V  ORDINATE  X  LO 

1  AO  OR  DISPLACEMENT  Y  LOAD  OR  DISPLACEMENT  PORE  PRESSURE  ) 

2004  FORMAT  !  I 1 2 , F 1 2 . 2 , 2F1 2 . 3 , 2E 24. 7, F 1 5. 3 » 

2005  FORMAT  ( 26H0N0DAL  POINT  CARD  ERROR  N#  15) 

2006  FORMAT (96H l ELEMENT  NO.  I  J  K  L  MATERIAL  X-ST 

♦RESS  Y-STRESS  XY-STRESS  Z-STRESS) 

2007  FORMAT! 112, 416, 1112, 4F12. 3) 
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2008  FORMAT  ( 29H0PRE SSURE  BOUNDARY  CONDITIONS/  40H  1  J 

♦  URE  I  PRESSURE  J 

2009  FORMAT  (2I6,2F12.3» 

2010  FORMAT  (12H1N.P.  NUMBER  IBX  2HUX  18X  2HUY  /  < 1 1 12, 2E20.7  ) 
END 


PRESS 

» 


r>  o 


j  — 

SUBROUTINE  STIFF 

l  IMPLICIT  REAL*8( A-HyO-Z I 

COMMON  ACELR* ACELZ * VOL f TEMP*  S I G I  I  400, 7) , HEDt 18) ,EI4. 121 »ROI 121 » 

I*R I  500 ) *  Z ( 500 ) , UR  1 500) ,UZ I  500 ) .CODE ( 500 ) y T 1 500 »  y PR! 100 y 2  )  , 

♦NUMMAT  t  NUMPC ,N ,MTVPE,KKK ,NUMNP  yNUMEL ,NNNy 
*  IBC  ( 100  )  y  JBC(  1001  yMTAGUOOl 

COMMON/ ARG/  RR I  5) ,ZZ I  5 1  y SI  10 , 10 1 y PI  8 ) y ST  I 3y 10 ) y C I  3, 3 ) , S  I  GIT) ,EE I  4 1 

I*,BOC lOOO)ySRlySR2« 

♦  RATIO! 400  I  ,LMI4),IXI400,5) yXCyYC 
COMMON  /BANARG/  MBANOyNUMBLK , B I  108 ) , AI 108 ,54 ) 


REWIND  2 
NB=27 
ND=2*NB 
ND2=2*ND 
ST0P=0.0 
NUMBLK*0 
DO  50  N= l ,N02 
B( N ) =0.0 
DO  50  M* l ,ND 
50  A  ( N ,  M )  =0  .0 

60  NUMBLK=NUMBLK+l 
NH=NB*I NUMBLKMI 
NM=NH-NB 
NL*NM-NB*1 
KSHI FT=2*NL-2 
DO  210  N= 1 »  NUMEL 
IF  ( I  X ( N , 5 ) )  210,210.65 
65  DO  80  1*1,4 

IF  (IXlNyll-NU  80,70,70 
70  IF  (IX(N.I) -NM )  90,90,80 
80  CONTINUE 
GO  TO  210 
90  CALL  QUAD 

IX(N,5)=-IXIN,5) 

IF(VOL)  100,100,110 
100  WRI T E 1 6 , 2000  )  N 
STOP* 1.0 

110  MM*4 

IFI IX(N,3)-IX(N,‘ ))  130,120,130 
120  MM*3 

130  DO  140  1=1, MM 
140  LMI I »«2*IX(N, I)-2 
DO  200  I = l , MM 
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! 


00  700  K-l.2 
1 1 -t HI  1  I ♦K-KSHIFT 
KK-7*I-2»K 
RtllhBIIIMPlKKl 
on  200  J-i *MM 
00  200  1-1*2 
JJ-LHtJlH“IIM*K  SHIFT 
U-2*J-2H 
1 F ( JJ I  200*200*175 
175  I Ff NO- JJ I  160*155* 155 
160  WRITE  16*2001)  N 
STOP-1. 0 
GO  TO  210 

156  AtlI*JJ}-A<II*JJ)*SIKK«LU 
200  CONTINUE 
210  CONTINUC 

00  220  N-Nl.NN 
K-?*N-K$H|f T 
fltK)  -6f  KUUZtN) 

St K- 11-81 K* l UUP  INI 
UNI-O. 

UZ INI-O. 

220  URtNI-O. 

IF  tNUNPCI  225.510.225 
225  nn  500  l-l.NUNPC 
I ■ IRCtl I 
J-JRCIU 
DR-Zt I )-Z( J) 

OZ-Rt Jl-Rt II 

PP2-tPRtL*2lPPRtL*l) 1/6. 

PP l-PP2*PRtl* 1 1/6. 
PP7-PP7*PR(L*2)/6. 

I l"2*l-KSHlfT 
J J-7* J-RSHIFT 
IF  til  I  265*265*215 
215  IFtll-NOI  260*240*265 
240  6m-l}-6tll-ll*PPl*0R 
Bill) -Mil I*PP140Z 
265  IFtJJI  100*100*270 
270  IFtJJ-NDI  275*275*100 
275  fttJJ-ll-ftf JJ-1UPP2*0R 
Rt  J.'l-Bt  JJUPP2P0Z 
100  CONTINUE 

510  00  400  --NL.NH 

IF  tM-NUWNPI  115*115*400 
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315  U"UMM) 

N»e*<-i-KSHIFT 

IF  ICODE(M) )  390,400*316 

316  IF  (COOEIMt-l • I  317*370,317 

317  IF  ICODEINI-2.)  310*390*318 

318  IF  (C00EIMI-3.I  390*380*390 
370  CALL  MOOIFY! A,B*N02*MBAND*N*UI 

CO  TO  400 

380  CALL  MOOIFY! A, B,N02,MBAN0*N*U) 

390  U»Ul!M| 

N«NM 

CALL  MODIFY! A*B*N02*MBAN0*N*UI 
400  CONTINUE 

C 

WRITE  121  IBIN),IAIN*M),M't*MBANOt*N*l,ND) 

C 

DO  420  N>1*N0 

K«N+NO 

B!NI*B|K) 

B!K)«0*0 
DO  420  M«l*ND 
AIN,M)>A!K,M) 

420  A! K  *MI »0.0 
C 

IF  INM-NUMNPI  60*400*400 
480  CONTINUE 
ACELR-O* 

ACELZ-O* 

NUMPC-0 

C 

IF! STOP  I  490*500.490 
490  CALL  EXIT 
500  RETURN 
C 

2000  FORMAT  I26H0NECATIVE  AREA  ELEMENT  NO.  141 

2001  FORMAT  I29HOOAND  WIDTH  EXCEEDS  ALLOWABLE  141 
ENO 
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SUBROUTINE  QUAD 
C 

IMPLICIT  REAL *6 1 A-H»  O-Z 1 

COMMON  ACE LR * ACEL2 • VOL* TEMP*  SIG1 (400* 7) *H€0( IB) *C<4* 121 *ROI 12) • 
*R( 300) *Z ( 500) * UR I  BOO  I • U2 ( 400  1 *COOE 1 500) • T 1 900 1  * PR I  ICO *2 1 • 

•NUMMAT *NUMPC*N*MTVPE*RRR*NUMNP«NUMEL »NNN» 

♦ IBC 1 1001 *JBC( 100) *MTACI400) 

COMMON/ARG/  RR 1 5)  *21(91* SI  10*  10)  *PIB)*ST(3*10)*CI3*3)*SIGI7)*EEI4) 
•  *B0( 10C0 )* SRI* SR? • 

•RAT  101400) «L M( 4) *  I XI 400*9) *XC*VC 
COMMON  /BANARG/  MBAN0*NUMBLK*B( 108) *AI 109*94) 

oimension  um.vni 

C 

CALL  STRSTR 
00  DO  J-1,10 
00  120  1-1*3 
120  ST  1 1  *J)«0* 

00  DO  1*1*10 
DO  S(l*J)*0. 
on  140  1*1*4 
NPP* I XI N* I ) 

RR 1 1  l«R(NPP) 

140  22 <  I  •  -2CNPP » 

XC*( RPC  I  )*RR( 2) +RR( 3)+RR(4) )/4« 

VC  *12/1  l)*Z2(2)*22(3)*22(4))/4« 

RRI9I-XC 

22I5I-VC 

K*9 

J«l 

1*4 

LMOI-9 

NT*4 

IF(IX(N.3>-IX(N*4))  160*190*160 
190  NT-1 

LMI 3 ) *9 
1-1 
K-3 
J-? 

XC*(RR( 1 • ♦RR( 2) ♦RRI 3))/3* 

VC-(Z2(l)+ZZ(2)*ZZ(3))/3* 

RRI9I-RRI3) 

22(91-22(3) 

160  00  200  NN-l.NT 
LM( 1)-2*|-1 
LMC  2)-2*J-l 
Ul 1 I -ZZ( J)-ZZ(K) 

U(2)-22(K)-22( I) 

•J(  3 1  *22 (  1  1-221  Jl 
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VI 1I*RRIK»-RRI  J» 

VI2I-PRII  l-RRIXI 
VI3)*RRI JI-RRIII 

ARf A*|RR| J|*UI2I+RRI 11*011 > ♦RR  €51 •U 1 1)1/2* 

VOL-VOL*ARFA 
COHM-.25/ARFA 
XNT-NT 
COM-2. 0/XNT 
COM*COM*COMM 

00  180  1*1*3 
I t*LM(L) 

STI  l*fIt*STIl*It)«UlLI*COM 

STI2*ttH)"STI2*t  IU I  ♦¥<!.)  *COM 

STI3*tl)*STI3*tt»*V(L»»COM 

Sri3*tt*ll*STt3*ttH»*UlU»COM 

00  180  M*l*3 

JJ*LMIM| 

SIII*«JJI*SIII*JJI*lUlLI*C< 1*  1  |4UIN|»VU. l*CI 3*  31 4VIMI»VILI4CI 1  *  31*U 
llM»»Ull»RCtl*3»4VINI|4C0MM 

SI ll*JJ*l)*S< Il.JJM  IMUll )*Ct 1 *2  t*Vt  M)*VIL )*CI 3*3)*UIMI ♦Vll l*CI 2» 
13l*Vl Ml tUILIRC 11*31 RUIMII 4C0MM 

SI  lt+l*JJ+ll*SII IM  *  JJM I  ♦!  VI  l  l*C  1 2*2|4VIM|»UIL |4CI3*3I6UIM|4UI L  •• 
1CI2*3I*VIM|AVIU*CI2*3I*UIM||»C0NM 
SI JJ»1*  I I l*SI 1 1 • JJ«1 1 
180  CONTINUE 
l*J 

200  CONTINUE 

IF!  IXIN.3I-IXIN.4H  220*280*220 
220  00  240  1*1*2 
KK-10-1 
00  240  K-|*KK 
CC*SIKK*1*K)/S(KK*1,KK»1) 

00  230  J*l*3 

230  STI J*K|-STIJ,K)-CC4S1IJ*KK4l) 

DO  240  J«1*KK 

240  SIJ*Kl*SI  J*K|-CC*SIJ*KIU1I 
250  CONTINUE 

ll*C 

IFINNN.EO.l)  GO  TO  280 
11*4 

260  S I G 1 1!*-SIGHN*I!M!»TEMP 
SIGI2)*-SIGI(N» 1 1^21 ♦TEMP 
SIGI3l*-StGIIN*IIO» 

00  520  l«l*8 
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pin  *0*0 
on  mo  j  - 1  •  3 

510  PI1I«PUI*STIJ*IIPS1GIJI 
520  PI  1 1 *P( I )*VOL 

OX»VOl*AC.f  IR*R0(MTYPE I  /4« 
OY-VOl*4CEU*RO<NTYPEI/4. 
00  530  I •  I  * 4 
PI  2* I >*Pi 2* 1 1 pOY 
530  P(2*l-l>*P*  2*l-mOX 
RETURN 
ENO 
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SUBROUTINE  STRSTR 


IMPLICIT  REAL*8lA-H.O-l>  11400.11  ,MEDtlBI.EIA.t2»**0t 12»* 

•IBCI 1001.  JBCC 1001.  5JJ5!4s?»l0aOI.PI»»*STI  3. 101 .Cl  »•»,  #SlCITIt««^» 
COMMON/ ARC/  RR13I.1H*)  *SI  10. 1»» 

♦,eOllOOO)*SAl.S*2. 


1»  l  XI  N»  1 1 

J- IXlN.2 I 

K-IXIN.3I 

L«|XlNfAI 

MTVPE-IXIN.5I 

VOL-O. 
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60 


TO 


tfmp.itiiiatui^tui^tilii/a.o 
00  50  KK-I*A 
EEIKKI*EIKK*MTVPEI 
IFIMTACINII  60.60*70 
CC-66€2»/«l*-66€?M 

RB-EEIII/U.-kf*  2JM2I 

C0MM*B3/II«-CC6*2I 

Cl  1 • 1 1 "COMM 
CI1.2»«C0MM*CC 
Cl  1*31*0# 

CI2. ll-CI 1.21 
CI2.2l-CII.il 
CI2.3I-0. 

CI3.lt -0. 

Cl  3.21-0. 

Cl 3.3I-.5*C0MM*I i.-CCI 

CC-0TANIEEI6I/5T.296I 

BB-DS0*TI9.0M2«09CC*CCI 

efiai-cc/bb 

FEI 3I-3.AEEI 31/BB 
GO  TO  500 

CC-DTANIEE I Al /5T.296I 
BB-DS0RTI9.0M2.0*CC*CCI 
EEIAI-CC/BB 

CC<«2,.*litMEI2ll/l3«-6.*EEI2l» 
00-ISIGUN.II-SIGIIN.2M/2. 
BJ2-I00*00^SIGIIN.3I*92I/I U  3* 

BJ2-DS0RT IBJ2I 


•  I EEI Al**2 1 1 
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BJl-l.S«(SIGI(N*!l»SIG!<N'?l >-3.  *E€  U>  *BJ2 
OD-0JI/BJ2 

6BM.*9.*(EEUI**2I*CC 
CC«3.*FE (4) RCC-DO/3. 

DD*EE(4l-D0/6. 

Hl«.5*CC/(&8*BJ2l 

H2-D0*CC/B8-EE(2I*EE(3I/<8B*6J2*;  t.-2.*EE(2) I  I 
H3« .5/1 B6*BJ2*BJ2I 

BB»f€tll/(l»*EEt2l) 

r.U,l»-BB*f  l.-H2-2.*Hl*SIGI  <N«  i  )-»H3*(SIGI  (N«  11**2)  I 
CUt?)--RB*(H2»Hl*tSIGI(Nt  l I *S IGI f N. 2) I *H3*S 1GI ( Nt i l*SIGI < N.2I I 
Cf  1.31  —  BB*(H|*SIGI(N»3i4H3*SIGItNtl)*SlGI(Nv3)l 
C(2.ll«C(1.2) 

C(2.2)«6B*U.-H2-2.*Hl*SlGI<Nt2l-H3*fSIGHN.2l**2ll 
C( 2. 3) — B8*IMl*S!GI f  N* 3) ♦H3*SIGI (N.2)*SIGI(N.3II 
Cf  3* 1 1 ■ C C t • 3 1 
Cf 3.2>«Ct2.3t 

C I  3 »  3 1*88*1 .5-H3*f  SI GlfN* 31**21) 

500  RETURN 
ENO 
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SUBROUTINE  MODIFY ( A, B*NEQ* MBAND* N *U> 
IMPLICIT  REAL*ft( A-H»0-Z ) 

DIMENSION  AI 108*94) *B( 108) 
no  250  M-2  tMBAND 
K-N-MM 

IF ( K I  235,235,230 
230  B(K)-BtK)-AIK*M)*U 
AIK, MI-0.0 
235  K«NtM-l 

IF(NEO-K)  250*240*240 
2 40  Rm-BtK)-A(N*Mt«U 
AIN, MI-0.0 
250  CONTINUE 
AIN* 1)-I*0 
BINI-U 
RETURN 
END 
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SUBROUTINE  BANSOl 
IMPLICIT  REAL*8I A-H»0-2 I 

COMMON  /BANARC/  MM,NUMBLK,B<  1061. A( 108,54) 
C 

NN«54 
Nt-NNM 
NH«NN+NN 
REM1N0  1 
RFWINO  2 
NB-0 

GO  TO  ISO 

C 

100  NB*NB*1 

DO  125  N*1 *NN 
NM«NN+N 
BINI-BINMI 
BINM|«0.0 
00  125  M*1»MM 
A(N,Ml»A(NM,MI 
125  AINM,M| »0»0 
C 

IF  (NUMBLK-NBI  150,200,150 
150  REAO  121  IBINI,IA(N,NI,N«1 ,MM) «N*Nl,NH) 

IF  INBI  200,100*200 

C 

20C  DO  300  N*1,NN 

IF  I  AIN, 1 1 )  225,300,225 
225  B( N I *BI Nl /Al N, 1 I 
DO  275  L*2*MM 
IF  lAIN.UI  230,275,230 
230  C-AIN, LI/AIN, II 
I*N*L-1 
J«0 

DO  250  K-l.MM 
J*JM 

250  Al |,J|-A||,JI-C*AIN,KI 
Bill »BI I l~AIN*L I *BINI 
A|N,U»C 
275  CONTINUE 
300  CONTINUE 
C 

IF  INUMBU-NBI  375,400,375 
375  WRITE  111  IBINI,IAIN,MI,M*2,MMI,N»1,NNI 
GO  TO  100 
C 

400  DO  450  M*1 ,NN 
N-NNM-M 
DO  425  K>2,MM 
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1-N4K-1 

4?5  MN»«B(N|-AIN,K)*B(U 
NM-NfNN 
BINM)«E<N) 

450  AINM.NB»*B(NI 
NB-NH-l 

IF  (NB)  475*500*475 
475  BACKSPACE  t 

READ  (1)  (R(N)*(A(N*M)*M»2*MM)*NM*NN) 
BACKSPACE  t 
GO  TO  400 

500  K*0 

UO  600  NB-UNUHBLK 
00  600  NM.NN 
K»K  +  l 
NM-N+NN 

600  R(K|*A(NM*NB) 

RETURN 

FNO 
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SUHMOUT INC  SIRESS 

IMPl  1C  I T  R£AL*0IA-H,O-2I 

COMMON  ACf LK,ACEL2,  VOL,  TfHP,  SIGH-00,?!  ,MCOI  IS)  ,H4*  121  ,«OI  12), 

•  M (5001, 2(500 1 ,UR<500 1, U2( 6001 «COOf 1 5001 ,113001 «PRl  ICO, 21, 

•NUMMAT  ,N(JMPC«N,MfVPC  ,KKK  ,NllMNP»NtJNFl,  ,NNN, 

•  I BC I IOO! *  JBCt ICO I  ,NT  AG (400 1 

COMMON/ ARC/  RRI5) *22151, SI  10, 101 .PI 01 ,ST| 3, 101 ,C(3,3I .SIGI 71 ,ECI4I 
• »BOI 1000) , SB I , SB?, 

•BAT  10(400 1  ,LMI4I,|X(400,5’,XC,VC 
COMMON  /BANARG/  MBANO,NUMBLK,BI 10B) ,A< 100,54 ) 

TOL-O. 

SR-1, 

NUMR-0 

MPRINT-0 

KK-0 

00  300  N-l.NUMEL 
RAT IGIN I - 1 • 

I XIN, 5 1 ■ I ABSI I X(N*5I  ! 

MTVPE-|X(N»5t 
CALL  QUAD 
MM-4 

IM  IXCN,  3t~lX(N,4tt  170,160,170 
160  MM-3 

170  00  1B0  1-1,3 
MB  1 1  1-0, 

00  100  J-l.MM 
II-26J 

JJ-261X1N, J1 

100  RR(II-PR(It*STII*!l!6B(JJt  +  ST(!,l!<”|)6B(JJ~l) 

00  190  1-1,3 

sicm-o. 

00  190  J-1,3 

190  SIG(ll-SIGm*CIl,Jt«RRUI 
00  195  1-1,3 
II-I64 

195  SIGI (N, 1 1 1-SIOIIN, 1 1  +  SIGI I ) 

00-1 SIGI (N,5t-SIGIIN,6I 1/2, 
AJ2«(00*006SIG1(N,7|*6P|/I1.-3.6|EE!4)662I) 

AJ2-0S0BTIAJ2I 

A Jl- 1,5*1 SIGI (N, 5 l»S I GIIN, 6 1 l-3.*EEI4l*AJ2 

r  AIL-AJ2+EE I4I*AJ1 

IFIMTAGINl.EO.O!  GO  TO  200 

|F(MTAG(Nt.FQ.2t  GO  TO  300 

OO-OABSI  f  A  IL -EE  I  31 1 

CHECK-. 02*EE13! 

IF ( OO-CHFCK 1  300,300,196 

196  KKK-1 
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(H'CHrcK/no 

IMfH-SW)  197,  100*  100 
197  SM*CH 

on  vn  300 

700  CONTINUE 

IFn-AlL.LT. EH3)  I  CO  TO  3C0 
KK-KK+I 

nO-tSIGI<N,l)-SIGI<N,2l)/2. 

BJ2-(00*DD»SIGHN,3I**2I/I l."3.*IFEI4l**2ll 
BJ2-DSQKTIBJ2I 

RJl-1.3*tSIGItN,l»*SIGI<N,?l  I-I.*EE  I4I*BJ2 

AAA-BJ2+BJ2 

BBft-AJ2*AJ2 

FFF«|.-3.*EE(4l*FE(4> 

CCO  tSIGI(N,l)-SIGI(N,2ll*(SIGI(N,3)-SIGI(N,6)l/4.*$IGI(N«3)*SIGl 
•IN, 71 

ODD-SIGI IN, 1 1 ♦SIGI (N.2I 
GGGaSIGI IN,3)*SIGI(N,6t 

ffmffffff 

GG-GGG*GGG 

00*000*000 

EF-2.23*EE(4)*EEI4) 

AA-AAA*FF-EF*D0 

BR-BBB*FF-EF*GG 

COCCC*FFF-EF*OOD*GGG 

00-1. 5*6E 141 •€! 131*000 

FF-1.5*EE(4t*EE<3l*GGG 

GG-EE(3I*EE<3) 

AAA-AA*B0-2.*CC 
RBB-AA-CC*nO-FF 
CCC-2.*00-GG*AA 
GGG-BBB*BBB-AAA*CCC 
GGG-OSORTIGGGI 
IFf AAAI  220,210,220 
210  RAT I0IN1 -• 5*CCC/BBB 
GO  TO  300 
220  AA»BBB/AAA 

BB-OABSIGGG/AAAI 

RATIOINI-AA-BB 

IF IRATIOINI.LT. 0.1  RAT|OIN)-AA*M 
IFIRATIOINI.GE.I.I  RAT IOINI -.99999 
IFIRATIOINI.LT. 0.1  RAT I0INI-0. 

300  CONTINUE 
C 

IFf KK.EQ.OI  GO  TO  420 
00  330  N-l,NUHEL 
IFIMTAGINI.GT.OI  GO  TO  330 
OO-RATIOINI 
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IF ( OO-SR I  305  *150?  350 

305  SR -00 
NUMR-N 
KKK*  I 

350  CONI INUE 

1FISR.LT. O.LI  SR-0.1 
IF INUNR.EQ.OI  GO  TO  520 

260  on  370  N*l,NUMEL 

l F ( MT  AGI N) • GT • 0 )  GO  TO  370 
IFIRATIOINI.LE.SRt  GO  TO  355 
00-RAT IOINI-SR 
IFI00-.05I  355*  355 .370 

355  MTAGINI-l 

370  CONTINUE 

420  CONTINUE 

00  410  N* .  * NUMNP 
I l«2*N-i 

fli  i  n«fli  ii  i*sr 

410  R(  f  I  ♦  1 1  *R(  IIM)*SR 
00  600  N-I.NUMEL 
I  *  I  X  (  N*  i  I 
J= I  X I N»  ? I 
K* I X I N»  3 1 
l  *  I  X I  N*  4 1 
MTVPE-1 XIN. 51 
IFIK.EO.U  GO  TO  460 
XC-IRI  I  I  ♦R  C  JI*R|KURIlll/4, 
YC*I  21  II^ZI  JI»2lKI*2lin/4. 
GO  TO  4  70 

460  XC-IRI I  I ♦R IJI^RIK I  1/3. 

VC=I  21  IUZI  JUZIK  11/3. 

470  CONTINUE 

00  450  1*1,3 
1  1*144 

SIGI  ll-SIGI(N,m-SlGI(N,ll 
SIGIIN, I ll*SIG( I >*l l.-SRI 
sigiin,  m*-siGitN,n  i 
SIGI  n*SIG(  l»*SR*SIGIIN,ll 

450  SIGIIN, I l-SIGI I  I 

SIGI7I«EF(2I*(SIGI II ♦SIGI 2  I  I 
DO  465  |*l»4 

455  Efc ( I  I *E ( I ,MTVPE I 

IFINTAGINI.EO.O)  GO  TO  500 
CC-0TANIEEI4I/57.296I 
0B*OSORTI9.O412.O*CC4CC1 
FF  1 4 ) *CC/BB 
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nn»i  siGm-siGun/?. 

RJ2=(0D*n0tSIG(3)**2)/(  l.-3.*(EE(4)**2) ) 

RJ2»nSQRT(RJ2) 

M(.|  /h.s*l\1(,(  I  MSI  G  )  -  l.  *Ff  ( <*  )  *RJ2 

SOO  SIM  (  N|  4  t  *  S  IM  n 

CC=(SIG(mSIG(2)  1/2. 

RR=(  SIGI 1 l-SIGI 21 ) / 2 • 

CR  =  0SQRT(R8*^2«,SIG(3I**2) 

SIG(4)*CC+CR 
S IGI 5 )*CC-CR 
S IG( 6 ) *0.0 

IF  ( (BB.FQ.O.O).AND. (SIG! 3).EQ.0.0)I  GO  TO  510 
SIG(6)=28.698*DAT  AN2 ( S IG( 3 ) *BB I 
510  CONTINUE 

IF(MPRINT)  520*520*550 
520  WRITE! 6*2000 )  NNN 
MPR I NT*50 

550  MPRINT=MPRINT-1 

WRITE  (6*  2001)  N,XC.VC.  (SIGUI  .  I  •1*71  ,HTAG(N) 

IF ( MTAGI N) • E0.0)  GO  TO  560 
DD=(  SIGI(N,5)-SIGI(N,6l )/2. 

B  J2*( 0D*00tSIGI IN*  71**2) /( l.-3.*(EE(4)**2l) 

R J2*  DSQRT  ( BJ2 ) 

T0L=T0L+BJ2 
GO  TO  600 

560  SIG(7)=EE(2)*(SIGI(N,5)*SIGl!N*6) ) 

BJ2*  OSQRT ( ( (SIGI (N,5)-SIGI(N,6) )  **2M  SI  GI  ( N  ,6  )-SIG(  7 1 )  **2  +  (  SIG  I  7 ) 
*-$ IGI (N,5) )**2)/6.0^SIGI(N, 7)6*2) 

T0L=T0L*BJ2 
600  CONTINUE 

SR2*  ( SR 1*SR )  ♦  SR2 
SRI*  ll.O-SR)  *  SRI 
WRITE (6  *2002 )  TQL.SR  ,NUMR,KK,SR2 
I F ( TOL-l • )  660,660,650 
650  KKK* 1 

GO  TO  700 
660  KKK*0 
700  CONTINUE 
11=0 

DO  710  1=1, NUMEL 
710  IFIMTAGUI.GE.il  H  *  1 1  ♦! 

IF(KK,GE. NUMEL)  CALL  EXIT 
800  RETURN 

2000  FORMAT! 1H1/ 

*36H  STRESSES  AFTER  APPROXIMATION  NUMBER  14//// 

*7H  EL. NO.  7 X  1HX  7X  IHY  4X  8HX-STRESS  4X  8HY-STRESS  3X  9HXY-STRESS 
*  2X  10HMAX-STRESS  2X  10HM IN-STRESS  7H  ANGLE  4X  8H2- STRESS  5X  7HPL 


58 


♦AST  I C  ) 

2C0  1  FORMAT  < I  7 , 2F8. 2, 1P5E 12 .4 . OP  IF 7 . 2 » IPE 12 .4 , 1  12 1 
2 COP  FORMAT ( 39H0THE  UNBALANCED  LOAO  A!  THIS  STAGE  IS  El*. 5// 

♦47H  THF  RATIO  FOR  CORRECTION  OF  STORED  STRESSES  IS  F10.4// 

* 3 1 H  THE  NEXT  ELEMENT  YIELDING  IS  I*/ 

*9 1 H  AND  THE  TOTAL  NUMBER  OF  ELEMENTS  THAT  CAN  YIELD  WITH  THE  LINEA 
*P  ADDITION  OF  TOTAL  LOAO  IS  14/ 

*50H  LOAD  UP  TO  THIS  STAGE  AS  A  FRACTION  OF  TOTAL  IS  F20.5  1 

C 

END 
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3.4  Kxamplo 


Naghdi  (19G7)  solved  the  problem  of  an  elastic-perfcctly  plastic  wcd|(c  under 
uniform  loading  on  one  face  (Fig.  3.2).  Plane  strain  conditions  were  considered. 
The  wedge  material  was  assumed  to  yield  according  to  Von  Mlscs*  yield  criterion. 
This  type  of  material  is  a  special  ease  of  Mohr-Coulomo  material  having  the  angle 
of  internal  friction  <t>  =  0 . 

Figures  3.3  and  3.4  show  the  theoretical  and  computed  results  for  the  dis¬ 
tribution  of  radial  and  circumferential  stress  at  various  stages  of  loadings.  The 
angle  t/>  denotes  the  angle  upto  which  the  yielding  has  progressed  from  the  boun¬ 
daries.  Fig.  3.5  shows  the  radial  strain  distribution  at  various  stages. 

Generally  the  agreement  between  results  computed  by  the  method  outlined 
and  the  exact  analysis  were  found  to  be  good. 
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CALCULATION  OF  STRESS  RATIO 


FIG  3.1 
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W  Ib/in 


FINITE  ELEMENT  IDEALIZATION  FOR  ELASTIC- PLASTIC  WEDGE 


FIG  3.2 
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DISTRIBUTION  OF  RADIAL  STRESS  FOR 
WEDGE  AT  VARIOUS  LOADS 


FIG  3.3 
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0 

DISTRIBUTION  OF  CIRCUMFERENTIAL  STRESS 
FOR  WEDGE  AT  VARIOUS  LOADS 


FIG  3.4 
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—  Noghdi's  theoretical 
solution 


o  Finite  element  results 


CHAPTER  IV 


COMPUTER  PROGRAM  FOR  ANALYSIS  OF  JOINTED  ROCK 


Chapter  IV,  Computer  Program  for  AniilyntH  of  Jointed  itoek 


4. 1  Organization 

The  computer  program  described  here  is  based  on  the  theory  described  In 
Chapter  I  and  II.  The  rock  mass  Is  considered  as  a  linear  elastic  material  in  the 
direction  of  compressive  stresses  and  is  assumed  to  have  no  resistance  to  defor¬ 
mation  in  the  direction  of  principal  stresses.  The  program  corrects  the  discrep¬ 
ancy  in  the  method  presented  by  Zienkiewicz  et  al  (1968).  This  was  pointed  out 
towards  the  end  of  Chapter  I.  It  makes  allowance  (or  the  boundary  conditions, 
residual  stresses,  stresses  due  to  temperature  change,  and  varying  pressure 
boundaries.  This  program  also  usrs  the  quadrilateral  elements  and  generates 
stiffness  the  same  way  as  that  described  in  Chapter  III. 

The  principal  program  called  MAIN  controls  all  the  data  input  and  control 
information.  It  does  the  system  initialization,  prints  the  control  data,  geometrical 
and  material  properties.  MAIN  calls  the  subroutines  for  stiffness,  solution  of 
equations  and  stress  calculations. 

4.12.  Stiffness  Matrix 

Stiffness  matrix  for  the  analysis  is  computed  in  blocks  by  the  subroutine 
STIFF.  For  element  stiffness,  it  calls  QUAD  for  triangular  and  quadrilateral 
elements  and  ONED  for  bar  elements.  Direct  stiffness  technique  was  used  to  get 
the  total  stiffness.  Equations  arc  modified  for  displacement  boundary  conditions. 
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4.13.  Load  Matrix 


Load  matrix  for  the  analysis  is  computed  in  LOAD  subroutine.  The  varying 
pressure  boundary  is  taken  into  account.  The  load  matrix  is  modified  in  the  LOAD 
subroutine  for  each  iteration  performed.  This  accounts  for  the  nonlinearity  intro¬ 
duced  by  progressive  cracking  by  considering  the  change  in  element  stiffness  as  a 
pseudo  load. 

4. 14.  Calculations  of  Displacements 

After  the  stiffness  and  load  matrices  for  a  stage  have  been  computed,  the 
resulting  equations  are  solved  by  calling  subroutine  SYMBAN.  This  subroutine 
was  developed  to  take  the  advantage  of  the  fact  that  stiffness  matrix  remains  the 
same  throughout.  Gaussian  elimination  technique  is  used  to  solve  the  equations. 

The  elimination  is  done  once  for  all  and  the  reduced  matrix  stored  on  auxiliary 
units.  Solution  for  each  iteration  consists  of  back-substitution  only.  This  approach 
results  in  considerable  economy  in  machine  time. 

4.15.  Calculations  of  Stresses 

After  the  displacements  have  been  computed,  the  stresses  are  computed 
using  the  constitutive  law.  A  change  in  the  minor  principal  stress  corresponding 
to  Poisson's  effect  has  been  introduced  when  one  of  the  principal  stresses  was 
compressive  and  others  tensile. 


67 


4.2.  Input  Data  Preparation 

1.  Control  Card  (A6).  This  card  will  carry  the  characters  START  in  columns 
1-5.  This  will  start  the  processing  of  the  data  deck  which  consists  of  the 
following  set  of  cards. 

2.  Job  Title  (72H).  This  card  will  give  the  descriptive  identification  for  the  job. 

3.  Control  Information  (415,  3F10.2,  15,  E15.4) 


Information  Columns 

Total  number  of  nodal  points  1-5 

Total  number  of  elements  6-10 

Number  of  different  materials  11-15 

Number  of  pressure  boundary  cards  Ifi  -  20 

Body  Force  in  X-direction  21-30 

Body  Force  in  Y-direction  31-40 

Reference  (stress-free)  temperature  41  -  50 

Number  of  Iterations  51-55 

Tolerance  to  Convergence  56-70 


4.  Material  Property  Cards.  One  set  of  cards  must  be  provided  for  each  material. 
In  each  set: 

a.  First  card  (215,  F10.3,  15)  will  give  the  following  information: 

Material  identification  number  1-5 

Number  of  temperature  cards  (8  maximum)  6-10 

Mass  density  of  the  material  11  -  20 
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21  -  20 


Material  code  to  designate  materials 
which  cannot  take  tension 

code  =  1  for  materials  which  cannot  take  tension 
0  for  materials  which  can  take  tension. 

b.  Subsequent  cards,  one  for  each  temperature,  the  number  being 

defined  in  columns  6-10  of  the  first  card,  will  carry  the  following  infor¬ 


mation  (4F10.0): 

Information  Columns 

Temperature  1-10 

Elastic  modulus  11  -  20 

Poisson's  ratio  21-30 

Coefficient  of  thermal  expansion  31  -  40 


5.  Nodal  Point  Cards  (15,  F5.0,  5F10.0) 

One  card  for  each  nodal  point  with  the  following  information.* 


Nodal  point  number 

1  -  5 

Type  of  nodal  point 

6-10 

X-ordinate 

11  -  20 

Y-ordinate 

21  -  30 

XR 

31  -  40 

XZ 

41  -  50 

Temperature 

51  -  60 

If  the  number  in  columns  6-10  is 

Zero  XR  is  the  specified  X-load  and  XZ  is  the  specified  Y-load 
1  XR  is  the  specified  X-displacement  and  XZ  is  the  specified  Y-load 
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2  XR  is  the  specified  X-load  and  XZ  is  the  specified  Y-displacement 

3  XR  is  the  specified  X-displacement  and  XZ  is  the  specified  Y-displacement 
All  loads  are  considered  to  be  total  forces  acting  on  an  element  of  unit  thickness. 
Nodal  point  cards  must  be  in  numerical  sequence.  If  cards  are  omitted,  the 
omitted  nodal  points  are  generated  at  equal  intervals  along  a  straight  line  between 
the  defined  nodal  points.  The  necessary  temperatures  are  determined  by  linear 
interpolation.  The  type  of  the  nodal  point,  as  well  as  XR,  XZ,  are  set  equal  to 
zero. 

6.  Element  Material  Cards  (1215) 

These  cards  shall  carry  the  material  type  of  all  the  elements.  Each  card 
shall  have  material  types  for  12  elements  in  sequence.  The  material  type 
for  each  element  must  be  read  in  as  no  interpolation  has  been  provided  for. 

7.  Element  Cards  (515,  5X,  3  F10.0) 

One  card  for  each  element  will  provide  the  following  data. 


Information  Columns 

Number  of  element  1  -  5 

Nodal  point  I  6-10 

Nodal  point  J  11-15 

Nodal  point  K  16-20 

Nodal  point  L  21-25 

Initial  stresses: 

(i)  component  in  x-direction  31  -  40 

(ii)  component  in  y-direction  41  -  50 

(iii)  shearing  stress  on  x-y  planes  51-60 
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8. 


Nodal  points  I,J,K,L  are  corners  of  each  individual  element  in  a  counter¬ 
clockwise  order  for  a  right  handed  system  of  coordinates.  For  triangular 
elements  set  nodal  point  L  same  as  nodal  point  K.  The  element  cards 
must  be  in  the  numerical  sequence.  Any  cards  that  are  omitted  will  be 
automatically  generated  in  the  program  by  incrementing  each  of  the  I,  J, 

K,  and  L  nodal  points  by  one. 

Pressure  Boundary  Cards  (215,  F10.0) 

One  card  for  each  boundary  element  which  is  subjected  to  a  normal 
pressure  will  carry  the  following  information: 

Information  Columns 


Nodal  Point  I 
Nodal  Point  J 
Normal  Pressure  at  I 


1-5 
6-10 
11  -  20 
21  -  30 


As  shown  in  the  sketch,  the  boundary  element  must  be  on  the  left  as  one 
progresses  from  I  to  J.  Surface  tensile  force  is  input  as  a  negative  pressure. 


Output  Information: 

The  following  information  is  developed  and  printed  by  the  program: 

1.  Reprint  of  input  data 

2.  Nodal  point  displacements 

3.  Stresses  at  the  center  of  each  element 
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4.3  Program  Listing 

* *4*  *<>;.  **  +  +****♦♦***********»***  +  **  ******t*ttt*t**t*«t«t*  +  tt  ****** 

*  fwn  DIMENSIONAL  ANALYSIS  Of  A  NO  TENSION  SYSTEM  * 

>:■*>  s*  **********  ******* **************  *****************  *********** 

nouii  PRECISION  S,C,B,A,P,ST,SIG,U,  V,CC,dB,CR,  AREA, COMM, DU, OV,DL, 
l DX , D f , X L  « RR  « FOR «  T DL «  UK  *1)7.  ,  S I  GI .  VOL ,  COM,  WOROI  2 1  tWORDl ,  E  ,COSA  ,  S INA 
2  *C'J ,  Z  Z  ,  C  l »  XC  » YC,  PR  *ACELR  ,  ACELZ  ,  TfcMP,  T,  R,Q,RO,Z 
t.DR.IJ/.DT.ZX 

r  d«m  IN  ACFLR,ACFl  /,TFMP,Q,ROI  12  I , H ( 9C0  I  *  2(900)  *  T  ( 90G  )*  PR  1 200*21 « 
1N!J'«NP,N'JMFI  , NUMMAT ,NUMPC,M TYPE, ILL, N,HED I  18 1  ,NTCI  10 1  *  COOE  f  9001 . 
2MK( 10) , IBCIPCCI , JBCI200),NCHfcCK 
CDMM()N/SYSA«G/UR(900)  tUZIOOOl  ,SIGII900.6),CUI  1600 )« TOL, VOL 
1  *  F  (  P  ,  4 » 121 

rnvM  )N  /ARG/C I  3 «  3 ) « S I 10* 10 ) , S IG( b I , PI  8 ) , STI 3, 10) » RRI 5) , 2Z( 5 ) t 
IXf ,YC,F FI  3) ,LM|4I , 1X1800,5) 

COMMON  /HANAKG/  MBANO,NUMBLK ,811800 ) ,A( 108*54) ,KKK,JA 
DATA  WDK0/6HST ART  ,6HST0P  / 

OFF  I NF  FILE  1 1 50 , | 500 ,U,NBK ) , ? I  800 , 256 ,U , I D) 

CAU  f PRSFTI2CH,256,-l,l) 

I  r  Kf AU( 4, IC06)  WOkOl 

iriwOROi.eo.woRoun  go  to  30 
I  r  i  mi  mo  1  •  eq*  woroi  2 1 )  stop 

GO  TO  10 

to  RFAO  15, 1000  I  HFO, NUMNP.NUMEL, NUMMAT, NUMPC , ACELR , ACELZ »  Q  ,NP  , TOL 
WMJ  TF<G,2000)  HEO,NUMNP,NUMEL, NUMMAT, NUMPC, ACELR,ACELZ,0,NP, TOL 
40  DO  V)  M-l, NUMMAT 

Pf  AD  (‘.,10011  MTYPF.NTCIMTYPF  ),Mn|MTYPE),MTC(  MTVPE  ) 
wM  II  14,2001)  MTYPF,NTCIMTYPE),ROIMTVPE),MTCIMTVPE) 

ND«K=NTC|MTYPfc  I 

R 1  Ai;  14,10021  (  If  I  I,  J,MTYPF.),J«l,4),  I*l,NUMTC) 
fcRITF  (4,20021  I (FI I,J,MTYPE) ,J-l,4> ,I«l,NUMTC) 

50  CONTINUE 

WPITF  (6,20031 
L  =0 

R I L ) =0  # 

7 (L ) =0. 

T ( L ) =0 . 

6'  P  F  A!)  (5,  1003)  N,CODF  IN)  ,R(N)*Z(N)«UR(N)*UZ(N)*T(N) 

F'L  =L  ♦  1 
Z  X  =  N-L 

l)R=(R(NI-RIL)  l/ZX 
DZ=( Z(N)-ZIL) l/ZX 
0T»( TINI-TILI  )/ZX 
70  L  =  !♦  I 

IF(N-l)  ICO, 90, BO 
pc  cnDf  (n=o.c 
Ml.l*HIL-ll*DR 
/II  )*ZIL-ll*nZ 
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T  ( I .)  *111.-1  )*DT 
ii«(U=r.n 
u/f  U-r.o 

GO  TO  70 

90  WHITE  (6*2004)  ( K , CODE ( K ) , P ( K ) , Z ( K ) ,UR ( K ) ,U2 ( K ) » T( K) * K»NL * N) 
1 F ( NUMNP-N)  100*110,60 
KO  WRITE  (6,2005)  N 
CAI.L  EXIT 
110  CONTINUE 

RF AO (5, 10071  (  IX(N,5)  »N=l,NUMEU 
WRITE  (6,2006) 

N  =  0 

130  RE  AD(  5,  1004 )  M,  (  I  X(  M ,  I )  ,  I  *  1 , 4)  ,  (  S  IGI  (  M ,  I )  ,  I * 1 , 3) 

IFCM.EO.il  GO  TO  140 
Z  X=M-N 

DO  135  1*1,3 

135  SIG(I)*(SIGI(M»!)-SIGI(N»!))/2X 
140  N*N+1 

IF  (M-NI  170,170,150 
15G  I  X  (  N  ,  1 )  *  I X  (  N“  1  ,  1 )  ♦  1 
IX(N,?)*IX(N-l«2)*l 
TX(N,3)*IX( N“l ,3)41 
IX(N,4)*IX(N- 1,41+1 
DO  160  1*1,3 

160  SI 01 (N,I)*SIGl(N-l, I ) ♦ S I G ( I ) 

170  WRIT  E(6,200 / )  N, ( I X( N, I) , I  * l , 5) , ( S IG I ( N, I ) , I  * l , 3) 

IF  (M-NI  180,180,140 
18C  IF  (NUMEl-N)  19U,190,130 
190  CONTINUE 

IF  (NUMPC)  290,310,290 
290  WRITE  (6,2008) 

DO  300  1*1, NUMPC 

RE  AIM  5, 1 CC5 )  IBC(L),JRC(L),PR(L,l),PR(L,2) 

300  WP I T E (6 , 20C9 )  I BC ( L I , JBC ( l ) ,PK< L, 1 1 , PR ( L , 21 
310  CONTINUE 

DO  440  N= 1 , NUMNP 

NN=2*N 

C.U  ( NN- 1  )  *0  • 

440  f U ( NN ) *0 • 

J*0 

on  340  N* 1 , NUMEL 
DO  340  1*1,4 
00  325  1*1,4 

KK=IABS( IX(N, I )-IX(N,L)) 

IF  (KK-J)  325,325,320 
320  J *KK 
326  CONTINUE 
340  CONTINUE 


73 


Wf'.AN  )-2*J+2 
WRITF|S,?01?)  MHANO 
3S0  CAI  L  STIFF 
NCHECK* 1 
K  KK  =  1 

TALL  SV MB AN 
K  K  K  =  2 

nn  S90  ILL=1,NP 

CALL  LOAD 

CAI  L  SYMBAN 

nil  400  N  =  1  ,  NUMNP 

MN-2*N 

CUINN-1 ) =CUINN-1 )*B(NN-1 ) 

400  CUINN)=CU(NN)+BINN) 

WRITEIS,2013)  LLL 

WRITEI6,2010)  I N ,fl 1 2*N- 1 ) , B I 2*NI *CU( 2*N- 1 ) vCU( 2*NI «  N  *  ItNUMNPI 
CALL  STRESS 

IF  INCHECK. EO.O)  GO  TO  600 
SOO  CONTINUE 
GO  TO  700 

600  WRITE  IS, 20 11 >  LI  l 
700  CO  TO  10 
999  FORMAT! IS) 

10CC  FORMAT  I 10A4/4I5,3F1O.2» I5,E15.4, 151 
10C1  FORMAT  I  215, 1F1C.3,I5I 

1002  FORMAT  I5F10.3) 

1003  FORMAT  I  I  5 , F 5 . 1 , 5F 10 . 4 ) 

1004  FORMAT  I  5  I  5 , 5X , 3F 10 . 4 ) 

1005  FORMATI2I5,2F10.3) 

1C 06  FORMAT  I A6) 

1007  FORMAT  I  1215) 

20 CO  FORMAT  I 1H1  10A4/ 


1  3 OHO  NUMBER  OF  NOOAL  POINTS -  13  / 

2  30H0  NUMBER  OF  ELEMENTS -  13  / 

3  30H0  NUMBER  OF  OIFF.  MATERIALS---  13  / 

4  30H0  NUMBER  OF  PRESSURE  CARDS-—  13  / 

5  30 HO  X-ACCELERAT  inN - E  12.4/ 

6  3CM0  Y-ACCELERATION - F.12.4/ 

7  30 HO  REFERENCE  TEMPERATURE - E12.4/ 

R  3CH0  NO.  OF  APPROXIMATIONS - 15/ 

9  30H0  TOLERANCE  FOR  CONVERGENCE—  E12.4I 


2001  FORMAT  I  1 7H0MATER IAL  NUMBER*  13,  30H,  NUMBER  OF  TEMPERATURE  CARDS* 
1  13,  15H,  MASS  DENSITY*  E12,4»16H,  MATERIAL  CODE*  151 

2002  FORMAT  I15H0  TFMPER ATURE  10X  5HEC  9X  6HNU  10X  5HALPHA/ 

1IF15.2,3E15.5) ) 

2CC3  FORMAT  I1C0H1NOOAL  POINT  TYPF  X  ORDINATE  Y  ORDINATE  X  LI) 

1 A 0  OR  DISPLACEMENT  Y  LOAD  OR  DISPLACEMENT  TEMPERATURE  I 
?CC4  rORMAT  I  I 12,F12.2,2F12.5,2E24.7»F12,3) 
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2005  FORMAT  (  26H0NOIJAL  POINT  CARO  fckROR  N  =  151 

?.nO(,  F  OHM  AT  (99H1FLFMFNT  NO.  I  J  K  L  MATERIAL 

1SIGIXX  SIC.IYV  SIGIXY  > 

20C7  FORMAT  (II l 3 ,4 16 , 1 1  12  * 3F 12 . 3  I 

2008  FORMAT ( 29HCPRESSURE  80UNDARY  CONDI T I ONS/40H  I  J  PRESSUR 

IF  I  PRESSURE  J) 

2009  F0RMAT(2I6,2F14.3) 

2010  FORMAT! I2H0N.P. NUMBER  17X  3HDUX  17X  3HDUY  I8X  2HUX  18X  2HUY / 

1(1 112.4F20.7I I 

2C11  FORMAT t  38FI0  NUMBER  OF  CYCLES  TO  CONVERGENCE  =  151 

2012  FORMAT(28H  BAND  WIDTH  FOR  THIS  DATA  =  15) 

2C13  FORMAT ( 30H1  RESULTS  OF  ITERATION  NO.'  15///) 

FNU 
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SURROUT INF  STIFF 


DOUBLE  PRECISION  S»C»B»A,P,S1»SIG,U»V»CC»BB,CR»  AREA » COMM, DU »DV*Dl, 
1DX,0Y,XL  ,PR,F0R,TC1L  ,UR,UZ  , S I G I , VOL ,COM , E ,COS A, S INA 
2,CU,ZZ,FF,XC,YC,PR ,ACELR  ,  ACELZ.TEMP, T,  RrQ»RO,2 
COMMON  ACELR.ACELZ ,TEMP,Q,RO( 12 1 «R( 900 l ♦ Z ( 900 1 « T 1900  I , PR (200  *21* 
INUMNP,  NIJMEL  »NIJMMAT  »NUMPC,MTYPE»LLL,N*HED(  18)  *  NTCC10)  , CODE (900)  « 

2 M T C (  10)  ,  IBCI200) , JBCI200) , NCHECK 
C0MM0N/SYSARG/URI900)  ,UZ(900) , SIGH  900, 6)  ,CU(  18001  ,TOL.VOL 
1 ,E(8,4,  12) 

COMMON  /ARG/C(3,3> ,S( 10,10) ,SIG(6),P(8),ST(3,10),RR(5)*ZZ(5) , 

1  XC,YC,FE! 3) ,LM(4) , IX( 800,5) 

COMMON  / BANARG/  MBAND, NUMBLK , B( 1800 ) ♦ Al 108, 54) ,KKK, JA 
DEFINE  FILE  1 (  50, 1 500, U, N8K  ) , 21 800 , 256, U , I D) 

NB  =  27 

ND=2*MB 

N02=2*ND 

ST0P=0.0 

MUMBLK=0 

JA=ND2*(MBAND+1 1/1500+1 
NBK=  l 

DO  50  N  =  1 » ND2 
DO  50  M=l,ND 
50  A(N,MI=0.0 
60  NUMBLK=NUMRLK+l 
NH=NB*( NUMBLK+1 ) 

NM=NH-NB  _  _  ....  .  _ 

NL  =  NM-NB+ l 
K  SH I FT=2*NL-2 
DO  210  N*l,NUMEL 
IF  (  I X ( N  ,  5 ) )  210,210,65 
66  DO  AO  1=1,4 

IF  ( I  X ( N, I )-NL )  AO »  70, 70 
70  IF  (  I  X ( N , I )-NM)  90,90,80 
80  CONTINUE 
GO  TO  210 

90  IF(IX(N,3)-IX(N,2) )  95,85,95 
85  CALL  ONED 

I  X  (  N ,  5 )  =- 1  X (N , 5 ) 

MM  =  2 

GO  TO  130 
95  CALL  QUAD 

IX(N,5)=-IX(N,5) 

I  0  =  N 

WRITEI2M0)  (  (C(  I  IK,  JJK)  ,  JJK-1,3)  ,EE  (II K  )  ,  1 1  K»l ,  3 ) 

1  , ( (S(JJI»KKI)»KK 1=1,8), JJI=1»8)»((ST( IKK, JKK) *  JKK=1 , 8) , IKK* 1,3) 

2  , (RR( J  I  I  ) ,ZZ( JI I ) ,JI 1  =  1,4)  ,XC,YC, TEMP, VOL 
IF(VOL)  100,100,110 
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ICO  WR I T F (ft  »  2000  1  N 

sinp«i.o 

IP  *M-4 

I  F I  I  <<*,  3)-IX(N,'.»  1  130.t2C.130 
120  MM=  3 

130  on  140  I  8 1 » MM 
140  I'M  1  1  *2* I  X ( N» I  1-2 
on  2 DC  l8l,MM 
on  200  K-1,2 

I I  ®LM(  I  H-K-KSHIFT 
KK=2*l-2+K 

on  200  J-l.MM 
00  200  L  8 1  ♦  2 
JJ=LM(J)*L-I  m-K SHIFT 
l L=2* J-2+L 
IF(JJ)  200.200,175 
175  IHNU-JJ)  180,195, 195 
180  WRITE  (6,2001  1  N 
STOP  *1.0 
GO  TO  210 

195  A(  I  I , JJ ) =A( I  I  ,JJ)*S(KK,LL) 

200  CONTINUE 
210  CONTINUF 

OH  4  DC  M*NL,NH 
IF  ( M-NUMNP 1  315,315,400 
315  U*UMM) 

N=2|M- 1 -KSH1  FT 
1 F ( CODE ( w  1  1  39C . 400,316 
31ft  I  F  (C'lOE  (  M  )-l  ,  )  317,370,317 
317  1 F  ( C'.DOE  (  Ml -  £•  1  318,390,318 
31«  ir(CniDF(M)-3.  1  390,380,390 
370  CALL  MOO  1 FY ( A ,N02  »M8AND, N ) 

GO  TO  400 

38C  CALL  MOO  I TY  ( A .N02.PBANP. N1 
390  U8UZ ( HI 
N*N  + 1 

CALI  MPO I F Y ( A»N02 ,MBANO»N 1 
400  CONTINUE 

WR  1  T  F  (  1  *  NRK  )  (  (MN,M),M81,MHAN01  ,N»i,NOI 
NBK*N8K ♦ J A 
00  420  N* 1 , ND 
K=N»ND 

00  420  M* l , NO 
A(N,M)»A(K,M) 

420  A(KfM)=0.0 

1MNM.L  T.NUMNP1  GO  TO  60 
nn  480  N= 1 , numfl 
I  X ( N , 5 )  =  I ABS ( IX(N,5) I 
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<*P0  CONTINUE 

I F I S  TOP  1  49C,4C0,49C 
40..  (.AIL  F  *  I  T 
r.00  RETURN 

?fCC  FORMAT  ( 76H0NEGAT I VF  AREA  ELEMENT  NO.  141 
2C01  FORMAT  ( 29H0BANI)  WIDTH  EXCEEDS  ALLOWABLE  141 
FND 
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:<tJHM'1lj|  INF.  MNfl) 

r. 

OnUBLC  PRECISION  S,C,  (3,A,P,  ST,  S  IG,U,  V,CC,RB,CR,AREA,COMM,OU»DV,OI., 
inX,l)YtXl.RR.FORtini  ,UR,UZ  #SIGI  ,VOl  .COM, l  ,COSA,S!NA 
2,CU,ZZ, Cf,XC,YC,PR,ACElR ,ACELZ,TEMP,T,  R,Q,RO,Z 
COMMON  ACElR»ACEL7»TEMP,Q,ROC 12 ) , R 1 900 ) , Z  <900 ) , TC 900 ) ♦ PRC  200 ,2 )  ♦ 
lNUMNP,NUMEL,NUMMAT,NUMPC,MTYPE,LLL,N,HEO< 18) ,NTCC 10 ) ,C0DEC900) ♦ 

2 MIC l 10 ) , IRC ( 200 ), JBC  <  200 ) «NC HECK 
COMMON/SYSARG/URC 900 ) , UZ ( 900 ) « S I G I C 900 ,6 ) tCUI  1800) ,TOL, VOL 
l  *  E  (  8 , 4 , 12) 

COMMON  / ARG/C ( 3*  3 ) *  SC 10, 10 ) , SIGC 6 ) , PC  8 ) , ST ( 3, 10) ,RRC 5) ♦ ZZ ( 5 )  « 
1XC,YC»EF<3),LM<4) ,  IXC800.5) 

COMMON  /BANARG/  MBANO, NUMBLK, BC  1800 ) , AC  108 ,54 ) ,KKK , JA 
C 

no  100  l»l,B 
P( I ) =0.0 
no  ioc  j=i,8 
IOC  S(I,J)=0.0 

MTYPE=IX(N,5) 

I s  I  X  C  N,  1  ) 

J* I  X ( N, 2 ) 

OX»R  C  J ) — R C  I  ) 

OY=ZCJ)-Z<  I  ) 

XL=DSQR  T ( OX**2+OY**2 ) 

COSA«OX/Xl 

SINAaDY/Xt 

COMM^EC l,2»MTYPE)*ECl,4,MTYPF)/XL 
C 

SC l , 1 ) »COSA*COSA*COMM 
SC l,2)»CnSA*SINA*C0MM 
S(l,3)*-S(  1,1) 

SC 1 » 4 ) *-SC 1,2) 

SC2, 1  )*S< 1,2) 

S  C  2 , 2 ) =S INA*S INA*CQMM 
SC2,3)*-SC 1,2) 

SC  2,4)=-S( 2,2) 

S  C  3 » 1 )*S ( l *3 ) 

SC3,?)=S<2»3) 

SC  3, 3)*S< 1,1) 

SC3,4)=S<1,2) 

SC  4, 1 )*S( 1 ,4) 

SC4,2)=SC2,4) 

SC  4, 3)*SC 3,4) 

SC4,4)=SC2,2) 

C 

C 

RETURN 

C 
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SUBROUTINE  QUAD 

DOUBLE  PRECISION  S,C,B,A,P,ST,SIG,U,V,CC,BB,CR«  AREA  *  COMM*  OUtOV  tOL  * 
ll)X,UY,Xl  ,RR,FUR,  TOL  ,UR,U7  ,S!GI ,VOL  ,COM,  E  ,  COS  A,  $  INA 
?*C«J*/ Z  •  cr  ,XC,YC,PR,ACELR,ACELZ,TEMP,T,  R,Q,RO*Z 
3,U.V,<T,XS.RAT10,0EN,XNT 

COMMON  ACELR,ACELZ,TEMP,Q,RG( 1 ? I , R (900 > , Z C 900 ) , TC90C ) .  PRC  200 . 2  ) . 
INIJMNP,  NUMEL  ,NUMMAT,NUMPC,MTYPF,LLL,N,HEOC  18 1  ,NTC  C  10 )  .CODE!  900) « 
?MTC(  10) ♦ I BC ( 2C  0 1 » JBC ( 200 ) » NCHECK 
C  OMMON/SYSARG/UR  <  900 1  »UZ(900)  «  S I G I  (900(6)  ,CUI  1000)  , TOL,  VOL 
1,E(R,4,12) 

COMMON  / ARG/C (3,3)  ,S(10*IQ),SIG(6),P(8)*ST(3,10)»RR(5),ZZ(5)» 

IXC  ,YC,FE ( 3) ,LM!4) , IX<800,5) 

COMMON  /BANARG/  MR AND, NUMBLK , 0 11800 ) ♦ A C 108 , 54 ) , KKK , JA 

niMFNSION  UC  3),VC  3) 

I  =  I  X  (  N,  1  ) 

J  =  I  X ( N, 2 ) 

K=IX(N,3) 

L  =  I  X ( N»4 I 
MTYPE=IX(N,5) 

VOL'O. 

TFMP*( T( I l*T( J)»T(K)+T(l ) )/4.C 

RATIO=G.C 
NUMTC=NTC(MTYPE ) 

IF  (NUMTC.FO.il  GO  TO  100 
00  50  M  =  2 , NUMTC 

IF  (  E(Mf  UMTYPFI-TEMP)  50,60,60 

50  CONTINUE 

60  l)FN  =  E(  M,  1,MTYPEI-E(M-1,  ltMTYPE) 

IF(OEN)  70,80,70 

70  RATIO-! TFMP-F(M-1 , 1 »MTVPE) ) /OEN 

BO  DO  90  K*  =  1 , 3 

90  FF(KK)-E(M-l,KK+l,MTYPE)+RATI0*(ECM,KK*i,MTYPE)-E4H-ltKK*lfMTYPE)) 

GO  TO  no 
ICP  DO  105  KK  =  1 , 3 
105  EF(KK)=E( l,KK+l,MTYPF) 

110  CONTINUF 

CPMM*EE( 1 )/ll.-EE(2)**2) 
r  <  1  ,  l)=rf)MM 
C ( l,?l=COMM*Fr (2) 

C  (  l  ♦  3 )  =  0  . 

C(2,n=C!l,2) 

f(2,2)*C(l,l) 

C  (  2 , 3 )  =  0 . 

r  ( 3,  n*o. 

C  (  3 , 2  )  =  0  • 

C( 3,3|=.5*COMM*(l.-EEC2) ) 
no  l 30  J=1,10 
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on  12"  1=1,3 
120  S T ( I , J)=0. 

nn  no  1=1,10 

130  S( i • j ) =o • 

DO  140  1=1,4 
NPP=IX(N,I> 

PR( I  l=R (NPP) 

140  7 / f  I  I =7  C  NPP> 

I F ( IX(N,3)-IX(N,4I)  145,150,145 
145  XC= ( RR( 1)+RK(2)+RR(3)+RR(4) 1/4, 

YC  =  (ZZm«7ZC?KZZ(3|4ZZ(4n/4. 

WR(5)=XC 
7  Z  f  5 1  =  VC 
K  =  5 
J=1 
1=4 

LM( 31=9 
NT  =  4 

GO  TO  160 
150  NT  =  1 

I  M( 31=5 
1  =  1 
K  =  3 
J  =  2 

XC  =  (RR(  1)+RR(2UPR(  3)1/3. 

YC=(ZZ( 1>+ZZ( 21 ♦ZZl 3)1/3. 

RR<5)=RR(3I 
Z  Z  (  5 »  =  Z  Z  (  3 ) 

160  DO  200  NN= 1 , NT 
LM(11=2*I-1 
IM( 2)=2*J-1 
U( 1 )*ZZ( J)-ZZ(K) 

U(?)  =  7Z(K1-ZZ(  1  ) 

Uf 31 =7  71  I  I — Z Z (  J ) 

V  (  1  ) =RR ( K I -RR ( J ) 

V( 2)*RR( I l-RP(K) 

V ( 3 ) =RP ( J l-RR ( 1  I 

ARPA=(RR(  J)*U(2)*RRm*um«>RR<5l*U<  311/2. 
vol=vol+area 

COMM= .75/ AREA 
XNT«NT 
r.0M=2./XNT 
f.  0M=C0M*C0MM 
DO  1 RO  1=1,3 
1 1 *LMI L ) 

ST(l,m*ST(l,  m+UU)*CO* 
ST(2»II+l)*ST(2»I  I+1)+V(L) *COM 
ST(3,in  =  ST(3,in+V(L»4C0M 


81 


ST  (3*IIM)  =  STI3*  1 1*1  MUtU*COM 
nn  irc  *=i.3 

S  (  I  I  , J J 1 =S1  I  I » J  J  )  ♦  ( :J  l  .  )*C<l,l)*UIM)*VtL)*C(3,3l*VtMU*C0MM 
SUI  »JJM)*SU  I,  JJ  +  l  I  HIM  l  )*CI1.2I*V(H)+VIL  )*Ct 3 t 3 )*UI M I) *COMM 
S(  I  1  »It  JJM  l*SII  IU,  JJ  +  1  )MVIU*Ct  l*  1  l*VI Ml  4-UI L  )*C (  3»3 1 *UI  M I  I  *COHH 
S  <  J  J ♦ 1 • I  I ) -S ( 1  I  *  JJ-M I 
1R0  CONTINUF 
1=1 
J  =  JM 
CONTINUF 

lFtIXIN,3l-lXtN,41l  220#250,220 
220  nn  240  1=1,2 
K  K  * l C“ I 
nn  240  K=1,KK 
CC=S<*KMtKi/S<KKM»KK  +  ll 
nn  2  jo  J  =  i,3 

23C  ST( J,K)«ST( J»K)-CC*STC J.KK+ll 
on  240  J= 1 , KK 

240  S(  J»K)  =  S(JfK  )-CC*SI  J»KKUI 
250  CONTINUF 
RF  TURN 
PNO 


SUBROUTINE  STRESS 

r 

™ H liar  PRFCISIUN  S*C»B»A«P,ST,SIG»U»V,CC»B&,CR»AREA,CnMM,DU»DV»DL, 
II)X,I)Y,X|.  ,RR,FfiR,Tni  ,UR»UZ  ,  S  1 GI  ♦  VOL  ,C  OM ,  E  » COSA  ♦  S INA 
2  ,CIJ ,  /  7  ,FE,XC»YC»PK,ACELR,ACELZ»TEMP»T  ,  R,Q,RO,Z 
3»CC»RH»CK*SS*SC,S2»C2»EPS,OT 

r.nMMQM  ACELR.  ACELZ,TEMP,Q,RO( 12>  ,R(90C  1,7(900)  ,  T(9GGI  ,PR(2G0,2I  t 
lNUMNP,NUMEL,NUMMAT,NUMPC,MTYPE,LLL,N,HEO(  18  1  ,NTC<  10  I  ,CODE  (  900 1  ♦ 
2MTC( 10) , IRC (200) , JBCI200) ,NCHECK 
CQMMON/SYSARG/UR(900) , UZ ( 900  I , S I G I ( 900 ,6 » ,CU(1800) ,TOL,VOL 
1 ,E(R,4, 121 

COMMON  / ARG/C 13*31 t S 1 10 . 1C  1 , S IG( 6 ) , P ( 8 1 , ST( 3, 101 »RR( 51 , ZZ ( 5 1  , 

1XC ,YC,FE(  31 ,LM( 41 , 1X1800,51 
COMMON  /RANARG/  M9AN0, NUMBLK , B I  1800 1 » A (  1 08 , 54 ) ,KKK , JA 
OFF  I NE  FILE  1(50, 1500, U,NBK 1 ,2(800,256, U, ID> 

FOR  =  0.0 
MPR I  NT  =  0 

00  600  M=  1 , NUMEL 
I  0=  M 

F  I N  0 ( 2 1  10) 

n=m 

IX(N,5I*IARS( I  X ( N , 5 )  I 
MTYPF=IX(N,5) 

S I G I ( N , 4  I =0 . 

SI  GUN,  5  1*0. 

SIGI (N,61*0. 

IT ( I X(N, 3)-IX(N,21)  90,60,90 
60  I  =  I  X ( N, 1  ) 

J  =  I  X ( N»  2  I 

XC=(R(I»*R( Jll/7.0 
YC«(Z( IlFZ(JII/?.0 
DX»R ( J ) -R (  I) 

0Y  =  /(,|»-/(  I  ) 

X(.  *  0  SORT  (DX**2*DY**2) 

mi=h(?*j-n-R(2*i-n 
nv=F (2*J1-R(2*t ) 

m  =r  v*i)Y/XL*ou*nx/xL 

SIG( 1)*E( 1 ,4,MTYPE l*DL*E( 1 ,2,MTYPEI/XLtSIGI(N,l )*E(1,4,MTYPE1 

I F ( S IG( 1 1 .GT ,0. 1  GO  TO  ICO 

SIGI(N,n»SIG(l) 

GO  TO  5C0 

100  Sir, I  (N,4)«E(l,2»MTYPE)*Dl/XOSIGI(N,ll 
SIGI (N, 11*0. 

GO  TO  420 

90  RFA0I2M0)  ((C(I!K,JJKl,JJK*l,3l,EE(IIK|,lIK*l,3l 

1  ,  (  (  S  (  J  J  I ,  K  K  1  )  ,KKI*l,8),JJI*l,8)  , (<ST(IKK,JKKI,JKK*l,8t,  IKK*  1,3  I 

2  , ( RR ( J I  I ) , ZZ ( JIM ,JI 1*1,41 ,XC,YC, TEMP, VOL 

MM=4 
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IFUXtM,31-IX(N,411  170, 1  60, 170 

L6C  =  3 

17C  on  me  i»i.3 
pp< i )=c. 
nn  me  j*i,mm 
i  i  =  2*j 

JJ=?*l X ( N,  J 1 

180  RR|  I )=RR(  I 1*ST(  I*  I  II*BI  I,  I  JJ-ll 

IF(LLL.GT.l)  GO  TO  102 
DT  =  TFMP-Q 
DX  =  E  £ ( 3  1 *0T 
DY  =  DX 

SIG<1)=-C(  l,ll*DX-C(l,2)*OV  ♦  S I G I  I N .  1 1 
SIG<2)*-C(2,l»*DX-CI2.2)*OY  ♦SIGI(N,2> 

S I G ( 3  1  *  S  I  G I (  N ,  3  ) 
on  TO  184 
18?  on  183  1*1,3 

183  S  l  G  ( 1 )=0.C 

184  CONTINUE 

00  190  1*1,3 
nn  i 8  s  j  =  i ,  3 

185  SIG(  I)*SIG(n»C(I,J)*RRlJ» 

190  CONTINUE 

IMIU.E0.1I  GO  TO  195 
DO  19?  1=1,3 

19?  S  I G I  I)*  S I  G(  I  1  +  S  I  G I  (  N ,  I  1 
195  CONTINUE 

CC  =  (  SIG{  l)*SIGI2n/2.0 
i3fi*l  SIC(  1  1-SIGI2) )/?. 

CP  *0S0R  T (80**2>S IG(3)**2) 

SIG(4)=CC4CR 
S I G ( 5  > -CC-CK 
S IG<  6  > -C  .0 

IE( (HH.FQ.O. ).ANn.(SIG(3).E0.0. ) )  GO  TO  200 
SIG(61=?8.648*CATAN2ISIG<3),BB} 

DX  =  ('  .0 

?  00  Sir, I  (N,  1  >  -  S  I G  4  l) 

SIGI (N,2)=SIG(21 
SIGI (N,3)=SIG(3I 

IF((SIG(4).LE.O.CO).OR.(MTC(MTVPE».EO.O»»  GO  TO  500 
IFISIGI5 1 .GE.O.CCOOl 1  GO  TO  370 
FPS=SIG (61/57,296 
CC  =  ['CnS(FPS) 

S  S  =  U  S I N ( EPS! 

C2=CC*CC 
S2*  SS*SS 

sc=ss*cc 
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PX  =  EF(2)*SIG(41 
S I G I <  M  *  A  >  =  SIG(4)*C2*DX*S2 
S1G I  IN*5»*  SIG(4)*S2*f3X*r2 
S I G I  (  N »  4 )  =  $  l  G  (  4  )  *  SC  -  D  X  *  SC 

r. 

r. 

GO  TO  400 

^70  SI(>I(N«4)*SIG(1) 

SIGI ( N. 5 ) =S IG( ?  ) 

SIGI (N»6) =SIG( 3) 

40C  SIGI(Ntn  =  SIG(l)"SIGI(N»4) 

SIGI (Nt2)»SIG(21 *SIGI (N,5> 

SIGI (Nt3)sSIG( 3 1  —  S I G I  IN, ft  I 
420  DX  =  S IG I (Nf4l**2*SIGI (N»5I**2+$IGI <N.6)**2 
DX=OSQRT(DXI 
IF(DX.LF.F0R1  GO  TO  450 
I  JK  =  N 
FOR  =  OX 

450  CONTINUE 

500  I F ( MPR I  NT )  55C, 520,550 
520  WR I T  E ( 6  t 2000  ) 

MPR I  NT  =  50 

550  MPRINT=MPPINT-1 

MP I T E (6  « 200 1IN«XC«VC*(SIG( I ) » I  *  1 »  6 ) »  DX 
400  CONTINUE 

WR I TE ( 6 » 2002 1  FOR  * IJK 
IFIFOR.LE.TDU  NCHECK  *  0 
RETURN 

20 CO  FORMAT  (7H1EL.N0.  7X  1HX  7X  IHV  4X  8HX-STRESS  4X  8HY-STRESS  3X 
1  9HXY-STPESS  2X  10HMAX-STRESS  2X  10HM IN-STRESS  7H  ANGLE  2X  17HUNB 
?AL  ANCEn  FORCF  1 

20C1  rnQM^r  (I7.2F8.2, 1P5F 12.4 , OP  IF  7. 2 f 1PE20.4I 

2002  FORMAT ( TPMGMAX I  MUM  UNBALANCED  FORCE  *  E12.5.16H  |N  ELEMENT  NO, 

t  15)  ’  oo 

F  NO 


su'M  out  iNf-  N*nn i r v (  a,nfq.mhano,n> 

IMillltLI  I'KK.ISUIN  A 
DIMENSION  A ( 108,541 
on  250  M  =  ? , MBANO 
K=N-M*1 

I F ( K . LF • 0  I  GO  TO  235 
A(K,M»=C.O 
?38  K=N+M-1 

IFtNEO.lT.K)  GO  TO  250 
A( N,M)=0.0 
250  CONTINUE 
A(N, 11=1 .0 
PFTURN 
f  ND 
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SURROUTINF  LOAO 

noURLE  PPFCIS1CN  S,C,B,A,P,ST,SIG,U,V,CC,i3B,CR, AREA.COMM.DU.DV.DL. 
lOX.DY.XI,  ,RP,FOR,TOlfUR,U?  .SIGI.VOL.COM.E.COSA.SINA 
?»CU,Z/,rF,XC»YC»PR»ACEtR»ACELZ»TEMP»T»  R.Q.RO.Z 
3.PP1 ,PP2,OR,D7,fP,XMM,DT 

COMMON  AC  FI  P , ACELZ , T EMP,  Q,RO( 12  I »RI 900 ) , Z 1900) »  T 1 90GI »PR I  200 .2 ) » 
l  NUMNP » NIIMFl ,  NUMMAT  »NUMPC  .HTYPE  iLll ,  N,  HEDC 18  I  ,NTCC  10 1  *  CODE  1 900)  » 
2MTCI 10 ) , I BC I  200) , J8C I  200 ) « NCHECK 
C0MM0N/SYSARG/URI900) ,UZI900I,S!GII900,6) ,CUC 1800) .TOL.VOL 
l.EC 8*4.121 

COMMON  /ARC/C  I  3. 3) , SI  10.10) • S1GI6 ) . P( 8 ) , ST  I  3 , 10) ,RRC5)*ZZ<5)  , 

1XC, YC.EF I  3) ,LMI4) ,  1X1800,5) 

COMMON  /RANARG/  MBANO, NUMBLK .BI1800 ) , A 1 108 , 54 1 ,KKK  ,  JA 
DO  50  N* 1 , NUMNP 
B ( 2*N- 1 ) =UR I N) 

R<2*NI=UZIN) 

UR  I N ) =0 • 

U7 I N ) -0 • 

50  CONTINUE 

IFIINUMPC.EO.O.OR.ILLL.GT.l)  )  GO  TO  300 
00  200  L  *  1 ,  NUMPC 
I  -  I HC I L  ) 

J  =  JRC I L  ) 

0R  =  7 I  I )-ZI J) 

D7  =  P I J ) -R I  I  ) 

PP2=|PR(L,2)«-PRIL,1 )  )/6. 

PP1*PP2+PR(L. l)/6. 

PP2=PP2*PRIL,2)/6. 

11=2*1 
J J=2* J 

HI  II-1)=RI I  1-1 )*PP1*0R 
B I II)=B1 II)*PP1*0Z 
Bl JJ-1)=RI JJ-1I*PP2*DR 
B(JJ)— BC JJ)+PP2*D7 
200  CONTINUE 
3C0  no  7C0  N=1 . NUMEL 
1*1X1 M.i) 

J  =  I  X I  N,  2 ) 

K  =  I  X  I  N»  3  ) 

L  =  I  X  |  N » 4 ) 

MTYPF=IXIN,5) 
rF(LU.EQ.lt  GO  TO  330 
IFIMTC(MTYPE) .EO.O)  GO  TO  700 
IFISIGI IN, 4) .NE.O. )  GO  TO  320 
IF(SIGI1N,5).NE.C.)  GO  TO  32C 
IF(StGIIN«6).NE.C.)  GO  TO  320 
GO  TO  700 
320  CONTINUE 
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Vi<"  IMJ.EP.KI  GO  TO  500 
II)  =  N 

P.r  AO(2‘  10)  (  <  C  (  I  IKtJJK)  « J JK=1 , 3)  , EH  I  IK)  •  I  IK* I,  31 

1  , ( <  S( JJI ,KKl»  ,KKI=l ,B> , JJI =1,8)  * (( ST ( IKK»JKK)»  JKK*l *  8 ) *  IKK* it  31 

2  • (RK(JI I ),Z7(JII),JI I *1, 4 1, XCt VC, TEMP, VOL 
IFCLLL.EQ.il  GO  TO  400 

S  I G  (  l)=-$IGHN,4> 

SIG( 2  I =-S I G I <N,5) 

SIG( 3)=-SIGI(N,6) 

GO  TO  450 
400  OT*TEMP-Q 
OX  =  f-F  (  3  )  *DT 
DY  =  E  fc ( 3  ) *0T 

S  I G  (  l)*-C(  1,1>*DX-Ct  1 » 2 ) *DY  +S  I G I  ( N»  1 ) 

SIG(2)*-CI2,1 )*OX-CI2,2)*OY  4-S1GI  (N,2) 

SIG(3)*SIGI(N,3I 
450  DO  520  1=1,8 
P( I >  =0 . 0 
1)0  510  J  =  1 »  3 

510  PUhPIll-STI  J,I)*SIGf  J) 

520  P(  I»=P(  I  I'WOL 

IF(lll.EO.l)  GO  TO  540 
00  530  1=1,3 
530  S IG( I >  =  0.C 
GO  TO  6C0 
540  mm=4 

IP(IX(N,3) ,EQ.IX(N,4))  HM=3 
XMM=MM 

r)Y=VOl*ACFl.Z*RO<MTYPE)/XMM 
nx=vm ♦AOELR*HOIMTYPE»/XMM 
no  550  1=1, MM 
P  (  2  *  I  >  =  P(2*mi)Y 
550  PC  2* I -1 1 «PC  2* I-l I ♦OX 
GO  TO  6C0 
500  CALL  0NED 

OX  =  R( J>-RI I  I 

PY=z(j)-zm 

F  P  =  - S l G l (  N ,  4 )  /  fc  ( l » 2 » MT  YPE 1 
I )  X  =  I)  X  *  fc  P 
0  Y  =  H  Y  *  f  P 

p(  n=su  ,i  )*nx«-$<  i , 2 1 *oy 

p(?)=5( 2,l)*OX+S(2,2)*OY 
P(3)=-P( 1) 

P  (  4  1  =  -  P  (  2  1 
GOO  no  6? 0  11=1,4 
520  |  M(  I  I  )  =  2*IX(N, lll-l 
no  650  J J=1 ,4 
II=IMUJ) 
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h(  1 1  i=ri  m+p(2*jj-n 
6*;o  nu mi»Bi 1 1 ♦  i )+p( 2*jj ) 

V'Q  CONTINUE 

nr  750  N=  1 » NUMNP 
I F  (COOL*  (  Nl  . E0.0.  I  GO  TO  750 
I F I ( CODE  INI .EQ.l. I .OR.(COOEIN) .EQ.3. 
If  I  I  CODE  INI . FQ. 2 . I .OR • ( CODE  I N) .EQ.3. 

750  CONTINUE 
RETURN 
FNO 


II  B  (  2*N- 1 1 «0» 
I)  BC2*N*«0.0 
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SUBROUTINE  SVMBAN 

untMLf  PRECISION  S,C,R,A,P,ST,S  IG»U» V»CC*BB»CR* AREA*COMM»DU»DV»DL t 
1DX,  r-V.XL  ,RR,FUR,Tnt,IJR,U2  »  S 1 G I  » VOL, COM, F  «COSA«SINA 
2  ,OJ ,  2Z , F£ , XC, YC,PR ,ACELR*  ACFL2 , TEMP, T,DT ,Q«RO 
r/lMHHN  /BANARG/  MB  ANO ,  NUMBLK  ,  B«  1  BOO  I , All  C6,  54 1  ,KKK  »  JA 
IJFFINF  FILF  1(  50 , 1 500  ,U  ♦  NRK I ,  2  (  BOO .  256  ♦  U  ♦  1 01 
C 

c 

NN*5S 

Nl*NNM 

NH*NN+NN 

NR=C 

K'NK*G 

NRK*  1 

FIND! I* 1  I 

IF(KKK.GT.l)  GO  TO  2000 
GO  TO  ISO 
IOC  NRK*NNK ♦ JA 
F INr ( l’NBKl 
Nfl*NB+ 1 
on  125  N*l,NN 
NM*NN*N 

00  125  M* 1 , MRAND 
A(N,M|«A(NM,M| 

125  A(NMtMI«0. 

r 

I F ( NOMRLK-NB )  150,200,150 
150  RFAOl  l’NBKHI  A(N,M|  ,  M«l,  MBANO)  »N*NL  »NH| 

NNK*NRK 

IFINRI  2C0,IC0,200 
200  DO  300  N*1 »NN 

|F(A(N,lll  225 »  3CC » 225 
225  00  275  1*2 ♦ MRAND 

I  F  (  A  (  N ,  L  I  I  230,275,230 
230  C*A(N, LI/AIN, 1) 

I=N+L-1 

J*0 

HO  250  K*L, MBANO 
J  =  JM 

250  A  ( I , J  I  =  A (  I»JI-C*A(N,K| 

A ( N, L I *C 
2  75  (.nNTINUE 
300  fONTlNUF 

NBK=NRK-JA 

WPITFC l »NRK  I  ( (AIN, Ml , M*1 , MBANO I ,N*1 ,NNI 
IFCNIIMRLK.FO.NBI  RETURN 

on  in  no 

200C  NQ*0 
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MFO=lHf f 
on  TO  4  5C 
4PC  NK  =  \fHl 

NHK*NNK4 JA 

FIND!  I'NBKI 
00  425  N-l »NN 
NM=NN*N 

PH  425  M=  l  f  BRAND 
A(N,M|=A(NM,M) 

425  A(NK,M)=0. 

IF(  MUMOLK-NBI  450,500,450 
450  P F  Al  I  I’NBKK  I A(N,M)  ,M»1 ,MBANDI , N=NL , NHI 
NNK*NBK 

I F I  MB »  50C  1 400  ,  500 
5C0  DO  55C  N* 1 » NN 
J*NQ+N 

DO  540  L*2,MttAND 
I = J+L- 1 

IFtNEQ-!l  545 , 540 ♦ 540 
540  HI  I  I *B( I »-A(N»LI*B( Jl 
545  IF(A(N,l)  .tO.O.J  AIN, 11  =  1. 

550  B( Jl =RI J)/A(N,1 I 
NHK=NNK-JA 

IF (NUBBLK.EQ.NBI  GO  TO  700 
600  NQ=NQ+NN 
GO  TO  4C0 
7^0  DO  750  M=  l , NN 
N=NNM-M 
J=MQ+N 

nn  750  L=2,MBAN0 
IF  (  A ( N«  L  I  I  740,750,740 
740  I=JH-1 

IF (NFO-I I  750,745,745 
745  FIJI  *H(  JI-AIN,LI*RU  I 
750  CONTINUE 
Mrt=NB- 1 

IFINft.EQ.C)  RF TURN 
FINOI I’NBK) 

DU  BOO  N*1 , NN 
NM*NN>N 

DO  800  M  =  1 , HBANO 
A(NM,M| = A ( N , M ) 

BOO  AIN,M)*0. 

READ  ll’NBKI  (I  A ( N, M I ,M* I , MBANDI , N= 1 ,NN) 

NBK=NBK-JA 

NO=NO-NN 

GO  TO  700 

FNO 
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4.4  Example 


A  tunnel  whose  configuration  and  finite  element  idealization  is  shown  in 
Fig.  4-1  was  solved  assuming  rock  to  crack  under  tension.  Lining  was  assumed 
to  be  capable  of  resisting  tensile  loads.  Figs.  4-2  and  4-3  show  the  initial  elastic 
and  final  'no  tension*  solutions.  The  tensile  stresses  are  indicated  by  arrows. 
Comparison  of  Fig.  4-2  and  4-3  shows  the  redistribution  of  stresses  caused  by 
the  inability  of  rock  to  withstand  tension.  In  the  example,  the  excavation  and 
lining  of  tunnel  is  assumed  to  be  a  single  step  ignoring  the  effect  of  sequential 
operations.  This  is  unrealistic  and  further  development  will  allow  for  actual 
sequence  of  construction.  However,  the  example  illustrates  use  of  the  computer 
program  and  is  similar  to  the  one  used  by  Zienkiewicz,  Valliappan  and  King  (1968). 

The  following  material  constants  were  used  for  the  solution: 

Lining:  E  =  2  x  10®  psi 

v  =0.16 

V  =  150  lbs. /ft. 3 

Rock:  E  =  1  x  10®  psi 

v  =  0.2 

Y  =  150  lbs. /ft. 3 
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FINITE  ELEMENT  IDEALIZATION  OF  LINED  TUNNEL 


ELASTIC  SOLUTION  FOR  THE  LINED  TUNNEL 
FIG,  4.2 
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CHAPTER  V 


ADDITIONAL  COMMENTS 


Chapter  V.  Additional  Comments 


The  work  reported  is  continuing  and  significant  changes  may  be  made 
before  the  finite  element  computer  programs  reach  their  final  form.  In  the 
case  of  elastic-plastic  materials  following  Mohr-Coulomb  Theory  non- monotonic 
loads  have  to  be  allowed  for  and  also  alternative  numerical  solution  procedures 
have  to  be  examined.  For  the  case  of  jointed  rock  the  program  included  in  the 
report  represents  a  necessary  first  step.  Modification  to  include  the  Griffith 
failure  criterion  is  being  done.  For  Griffith  rupture  caused  by  a  tensile  stress 
field  acting  on  pre-existing  flaws,  the  solution  scheme  appears  to  be  fairly 
straight  forward.  However,  the  case  of  fracture  under  compressive  stress 
fields  may  require  development  of  new  methods. 


List  of  References 


1.  Baker,  L.E.,  Sandhu,  R.S.  and  Shieh,  W.Y.  (1969),  "Application  of  Elasto- 

Plastic  Analysis  in  Rock  Mechanics  by  Fintie  Element  Method,"  Proc. 
Eleventh  Symposium  on  Rock  Mechanics,  Berkeley,  California 

2.  Bieniawski,  Z.T.  (1967),  "Mechanism  of  Brittle  Fracture  of  Rock,"  Parts 

I,  II,  in,  Int.  J.  Rock  Mech.  Min.  Sci.,  V.  4,  p.  395-525. 

3.  Brace,  W.  F. ,  (1964),  "Brittle  Fracture  of  Rocks,"  p.  m,  in  State  of  Stress  in 

the  Earth's  Crust,  ed.  W.  R.  Judd,  Am.  Elservier  Publishing  Co,  Inc. 

New  York. 

4.  Brace,  W.  F.  (1966),  "Dilatancy  in  the  Fracture  of  Crystalline  Rocks,"  J. 

Geoph.  Res.,  V.  71,  No.  16,  3939-53. 

5.  Brady,  B.T.  (1970),  "A  Mechanical  Equation  of  State  for  Brittle  Rock,  Part  I, 

The  Pre-Failure  Behavior  of  Brittle  Rock,"  Int.  J.  Rock  Mech.  Min.  Sci., 

V.  7,  385-421. 

6.  Brady,  B.T.  (1969a),  "The  Nonlinear  Mechanical  Behavior  of  Brittle  Rock, 

Part  I,  Stress-Strain  Behavior  During  Regions  I  and  n, "  Int.  J.  Rock  Mech. 
Min.  Sci.,  V.  6,  No.  2,  211-225. 

7.  Brady,  B.T.  (1969b),  "The  Nonlinear  Mechanical  Behavior  of  Brittle  Rock, 

Part  II,  Stress-Strain  Behavior  During  Regions  m  and  IV,  Int.  J.  Rock 
Mech.  Min.  Sci.,  V.  6,  No.  3,  301-310. 

8.  Clough,  R.W. ,  (1960),  "The  Finite  Element  Method  in  Plane  Stress  Analysis," 

Proc.  2nd  ASCE  Conf.  on  Electronic  Computation,  Pittsburgh,  Pa. 

9.  Clough,  R.  W.  (1965),  "The  Finite  Element  Method  in  Structural  Mechanics," 

Chapter  7,  STRESS  ANALYSIS,  ed.  O.C.  Zienkiewicz  andG.S.  Holister, 
London:  Wiley. 

10.  De  Arantes  e  Olivereira,  Eduardo  R.  (1968),  "Theoretical  Foundations  of  the 

Finite  Element  Method,"  Int.  J.  Solids  Structures,  V.  4,  929-952. 

11.  Drucker,  D.  C,  (1951),  "A  More  Fundamental  Approach  to  Plastic  Stress-Strain 

Relation,"  Proc.  First  U.S.  National  Congress  of  Applied  Mechanics,  ASME, 
487-491. 

12.  Drucker,  D.C.  and  Prager,  W.  (1952),  "Soil  Mechanics  and  Plastic  Analysis  or 

Limit  Design,"  Q.  App.  Math,  v.  10,  157-165. 


97 


13.  Drucker,  D.  C.,  Gibson,  R.  E.  and  Henkel,  D.  J.  (1955),  ’’Soil  Mechanics 

and  Work-Hardening  Theories  of  Plasticity.”  Paper  No.  2864,  Trans., 

ASCE. 

14.  Duncan,  J.M.  and  Goodman,  R.E.  (1968),  ’’Finite  Element  Analysis  of  Slopes 

in  Jointed  Rock,”  Report  No.  TE-6R-1,  U.S.  Army  Engineer  Waterways 
Exp.  Station,  Corps  of  Engineers. 

16.  Einstein,  H.H.  ,  Bruhn,  R.W.  and  Hirschfeld,  R.C.  (1970),  ’’Mechanics  of 

Jointed  Rock,  Experimental  and  Theoretical  Studies,”  Soil  Mech.  Publication 
No.  268  M.  I.  T. 

16.  Felippa,  C.  (1966),  ’’Refined  Finite  Element  Analysis  of  Linear  and  Nonlinear 

Two-Dimensional  Structures,  ”  Ph.  D.  Thesis,  University  of  California, 

Berkeley. 

17.  Felippa,  C.A.  and  Clough,  R.W.,  "The  Finite  Element  Method  in  Solid  Mechanics,” 

V.2,  SIAM-AMS  Proc.,  Numerical  Solution  of  Field  Problems  in  Continuum 
Physics,  Am.  Math.  Soc.,  Providence,  R.I, 

18.  Goodman,  R.E. ,  Taylor,  R.L.  and  Brekke,  T.L.  (1968),  "A  Model  for  the 

Mechanics  of  Jointed  Rock,”  ASCE,  SM3,  637-658. 

19.  Green,  A.E.  and  Naghdi,  P.M.  (1965),  "A  General  Theory  of  an  Elastic- 

Plastic  Continuum,”  Arch.  Rational  Mech.  Anal,  V.  18,  251-281. 

20.  Hill,  R.  (1950),  The  Mathematical  Theory  of  Plasticity.  Oxford  University 

Press. 

21.  Holand  and  Bell  (1970),  ’’Finite  Element  Methods  in  Stress  Analysis,”  TAPER, 

The  Technical  University  of  Norway,  Trondheim,  Norway,  2nd  print. 

22.  Koiter,  W.T.  (1953),  ’’Stress-Strain  Relations,  Uniqueness  and  Variational 

Theorems  for  Elastic- Plastic  Materials  with  a  Singular  Yield  Surface,” 

Quar.  Appl.  Math.  Vol.  11,  350-354. 

23.  Malina,  H.  (1970),  "The  Numerical  Determination  of  Stresses  and  Deformations 

in  Rock  Taking  into  Account  Discontinuities,”  Rock  Mechanics  2,  1-16. 

24.  Melosh,  R,  J.  (1963),  "Basis  for  Derivation  of  Matrices  for  the  Direct  Stiffness 

Method,”  AIAA,  V.  1,  No.  7. 


98 


25.  Naghdi,  P.M.  (1960),  "Stress -Strain  Relations  in  Plasticity  and  Thermo¬ 

plasticity,  "  Plasticity:  Proceedings  of  the  2nd  Symposium  on  Naval 
Structural  Mechanics,  Pergman  Press,  121-169. 

26.  Oden,  J.T.  (1969),  "A  General  Theory  of  Finite  Elements,  I.  Topological 

Considerations;  n.  Applications,  Int.  J.  Num.  Math.  Eng.,  V.  1, 

205-221,  247-259. 

27.  Plan,  T.H.H.,  (1968),  "Variational  Principles  and  Their  Application  to 

Finite  Element  Methods,"  Lecture  Notes  on  "Finite  Element  Methods 
in  Solid  Mechanics,"  24-38,  M.I.T. 

28.  Prager,  W.  (1949),  "Recent  Developments  in  the  Mathematical  Theory  of 

Plasticity,"  J.  App.  Phys.,  V.  20,  p.  235. 

29.  Proc.  Conf.  Matrix  Methods  in  Structural  Mechanics,  Wright- Patterson 

Air  Force  Base,  Ohio ,  1965. 

30.  Proc.  2nd  Conf.  on  Matrix  Methods  in  Structural  Mechanics,  AFFDL-TR- 

68-150,  Wright  Patterson  AFB,  Ohio,  1968. 

31.  Reyes,  S.F,  (1965),  "Elasto- Plastic  Analysis  of  Underground  Openings  by 

the  Finite  Element  Method,"  Ph.D.  Thesis,  University  of  Illinois. 

32.  Reyes,  S.F.  and  Deere,  D.  V.  (1966),  "Elastic- Plastic  Analysis  of  Under¬ 

ground  Openings  by  Finite  Element  Method,"  Proc. ,  First  Congress 
International  Society  of  Rock-Mechanics,  Lisbon. 

33.  Sandhu,  R.S.  and  Wilson,  E.L.  (1970),  "Finite  Element  Analysis  of  Stresses 

in  Mass  Concrete  Structures,"  Preprint,  Symposium  on  Impact  of  Computers 
on  the  Practice  of  Structural  Engineering  in  Concrete,  ACI,  St.  Louis, 
Missouri. 

34.  Swanson,  S.R.  (1970),  "Development  of  Constitutive  Equations  for  Rocks," 

Ph.  D.  Thesis,  University  of  Utah. 

35.  Walsh,  J.  B.  (1965a),  "The  Effect  of  Cracks  on  the  Compressibility  of  Rock, " 

J.  Geophys.  Res,,  V.  70,  No.  2,  381-389. 

36.  Walsh,  J.  B.  (1965b),  "The  Effect  of  Cracks  on  the  Uniaxial  Elastic  Com¬ 

pression  of  Rocks, "  J.  Geophys.  Res.,  V.  70,  No.  2,  399-411. 


99 


37.  Walsh,  J.  B.  (1965c),  "The  Effect  of  Cracks  in  Rocks  on  Poisson's  Ratio," 

J.  Geophysical  Res. ,  V.  70,  N.  20  ,  5249-57. 

38.  Wilson,  E.L.  (1963),  "Finite  Element  Method  in  Two-Dimensional  Problems," 

D.  Engg.  Thesis,  University  of  California,  Berkeley. 

39.  Wilson,  E.L.  (1965),  "Structural  Analysis  of  Axisymmetric  Solids,"  J.  A1AA, 

V.  3. 

40.  Yamada,  Y. ,  Yoshimura,  N.,  and  Sakural,  T.  (1968),  "Plastic  Stress-Strain 

Matrix  and  its  Application  for  Solution  of  Elastic- Plastic  Problems  by  the 
Finite  Element  Method,"  Int.  J.  Mech.  Sci.  ,  Pergamon  Press,  V.  10, 
343-354. 

41.  Zienkiewicz,  O.C.  and  Cheung,  Y.  K.,  (1967),  The  Finite  Element  Method 

in  Structural  and  Continuum  Mechanics.  McGraw-Hill. 

42.  Zienkiewicz,  O.C.,  Valliappan  and  King,  I. P.  (1969),  " El asto- Plastic 

Solutions  of  Engineering  Problems,  Initial  Stress,  Finite  Element  Approach," 
International  Journal  of  Numberical  Methods  in  Engg. ,  V.  1,  75-100. 

43.  Zienkiewicz,  O.C.,  Valliappan,  S. ,  and  King,  I.  P.  (1968),  "Stress  Analysis 

of  Rock  as  a  'No  Tension'  Material,"  Geotechnique  18,  No.  1,  56-66. 

44.  Zlamal,  M.  (1968),  "On  the  Finite  Element  Method,"  Numer.  Math.  12, 

394-409. 


100 


